| Geant4 Cross Reference |
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 // class G4ReplicaNavigation Implementation
26 //
27 // Author: P.Kent, 1996
28 // --------------------------------------------------------------------
29
30 #include "G4ReplicaNavigation.hh"
31
32 #include "G4AffineTransform.hh"
33 #include "G4SmartVoxelProxy.hh"
34 #include "G4SmartVoxelNode.hh"
35 #include "G4VSolid.hh"
36 #include "G4GeometryTolerance.hh"
37
38 namespace
39 {
40 const G4ThreeVector VecCartAxes[3]=
41 { G4ThreeVector(1.,0.,0.), G4ThreeVector(0.,1.,0.), G4ThreeVector(0.,0.,1.) };
42 const G4ExitNormal::ESide SideCartAxesPlus[3]=
43 { G4ExitNormal::kPX, G4ExitNormal::kPY, G4ExitNormal::kPZ };
44 const G4ExitNormal::ESide SideCartAxesMinus[3]=
45 { G4ExitNormal::kMX, G4ExitNormal::kMX, G4ExitNormal::kMX };
46 }
47
48 // ********************************************************************
49 // Constructor
50 // ********************************************************************
51 //
52 G4ReplicaNavigation::G4ReplicaNavigation()
53 {
54 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
55 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
56 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
57 halfkCarTolerance = kCarTolerance*0.5;
58 halfkRadTolerance = kRadTolerance*0.5;
59 halfkAngTolerance = kAngTolerance*0.5;
60 fMinStep = 0.05*kCarTolerance;
61 }
62
63 // ********************************************************************
64 // Inside
65 // ********************************************************************
66 //
67 EInside
68 G4ReplicaNavigation::Inside(const G4VPhysicalVolume* pVol,
69 const G4int replicaNo,
70 const G4ThreeVector& localPoint) const
71 {
72 EInside in = kOutside;
73
74 // Replication data
75 //
76 EAxis axis;
77 G4int nReplicas;
78 G4double width, offset;
79 G4bool consuming;
80
81 G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
82
83 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
84
85 switch (axis)
86 {
87 case kXAxis:
88 case kYAxis:
89 case kZAxis:
90 coord = std::fabs(localPoint(axis))-width*0.5;
91 if ( coord<=-halfkCarTolerance )
92 {
93 in = kInside;
94 }
95 else if ( coord<=halfkCarTolerance )
96 {
97 in = kSurface;
98 }
99 break;
100 case kPhi:
101 if ( (localPoint.y() != 0.0)||(localPoint.x() != 0.0) )
102 {
103 coord = std::fabs(std::atan2(localPoint.y(),localPoint.x()))-width*0.5;
104 if ( coord<=-halfkAngTolerance )
105 {
106 in = kInside;
107 }
108 else if ( coord<=halfkAngTolerance )
109 {
110 in = kSurface;
111 }
112 }
113 else
114 {
115 in = kSurface;
116 }
117 break;
118 case kRho:
119 rad2 = localPoint.perp2();
120 rmax = (replicaNo+1)*width+offset;
121 tolRMax2 = rmax-halfkRadTolerance;
122 tolRMax2 *= tolRMax2;
123 if ( rad2>tolRMax2 )
124 {
125 tolRMax2 = rmax+halfkRadTolerance;
126 tolRMax2 *= tolRMax2;
127 if ( rad2<=tolRMax2 )
128 {
129 in = kSurface;
130 }
131 }
132 else
133 {
134 // Known to be inside outer radius
135 //
136 if ( (replicaNo != 0)||(offset != 0.0) )
137 {
138 rmin = rmax-width;
139 tolRMin2 = rmin-halfkRadTolerance;
140 tolRMin2 *= tolRMin2;
141 if ( rad2>tolRMin2 )
142 {
143 tolRMin2 = rmin+halfkRadTolerance;
144 tolRMin2 *= tolRMin2;
145 if ( rad2>=tolRMin2 )
146 {
147 in = kInside;
148 }
149 else
150 {
151 in = kSurface;
152 }
153 }
154 }
155 else
156 {
157 in = kInside;
158 }
159 }
160 break;
161 default:
162 G4Exception("G4ReplicaNavigation::Inside()", "GeomNav0002",
163 FatalException, "Unknown axis!");
164 break;
165 }
166 return in;
167 }
168
169 // ********************************************************************
170 // DistanceToOut
171 // ********************************************************************
172 //
173 G4double
174 G4ReplicaNavigation::DistanceToOut(const G4VPhysicalVolume* pVol,
175 const G4int replicaNo,
176 const G4ThreeVector& localPoint) const
177 {
178 // Replication data
179 //
180 EAxis axis;
181 G4int nReplicas;
182 G4double width,offset;
183 G4bool consuming;
184
185 G4double safety = 0.;
186 G4double safe1,safe2;
187 G4double coord, rho, rmin, rmax;
188
189 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
190 switch(axis)
191 {
192 case kXAxis:
193 case kYAxis:
194 case kZAxis:
195 coord = localPoint(axis);
196 safe1 = width*0.5-coord;
197 safe2 = width*0.5+coord;
198 safety = (safe1<=safe2) ? safe1 : safe2;
199 break;
200 case kPhi:
201 if ( localPoint.y()<=0 )
202 {
203 safety = localPoint.x()*std::sin(width*0.5)
204 + localPoint.y()*std::cos(width*0.5);
205 }
206 else
207 {
208 safety = localPoint.x()*std::sin(width*0.5)
209 - localPoint.y()*std::cos(width*0.5);
210 }
211 break;
212 case kRho:
213 rho = localPoint.perp();
214 rmax = width*(replicaNo+1)+offset;
215 if ( (replicaNo != 0)||(offset != 0.0) )
216 {
217 rmin = rmax-width;
218 safe1 = rho-rmin;
219 safe2 = rmax-rho;
220 safety = (safe1<=safe2) ? safe1 : safe2;
221 }
222 else
223 {
224 safety = rmax-rho;
225 }
226 break;
227 default:
228 G4Exception("G4ReplicaNavigation::DistanceToOut()", "GeomNav0002",
229 FatalException, "Unknown axis!");
230 break;
231 }
232 return (safety >= halfkCarTolerance) ? safety : 0;
233 }
234
235 // ********************************************************************
236 // DistanceToOut
237 // ********************************************************************
238 //
239 G4double
240 G4ReplicaNavigation::DistanceToOut(const G4VPhysicalVolume* pVol,
241 const G4int replicaNo,
242 const G4ThreeVector& localPoint,
243 const G4ThreeVector& localDirection,
244 G4ExitNormal& arExitNormal ) const
245 {
246 // Replication data
247 //
248 EAxis axis;
249 G4int nReplicas;
250 G4double width, offset;
251 G4bool consuming;
252
253 G4double Dist=kInfinity;
254 G4double coord, Comp, lindist;
255 G4double signC = 0.0;
256 G4ExitNormal candidateNormal;
257
258 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
259 switch(axis)
260 {
261 case kXAxis:
262 case kYAxis:
263 case kZAxis:
264 coord = localPoint(axis);
265 Comp = localDirection(axis);
266 if ( Comp>0 )
267 {
268 lindist = width*0.5-coord;
269 Dist = (lindist>0) ? lindist/Comp : 0;
270 signC= 1.0;
271 }
272 else if ( Comp<0 )
273 {
274 lindist = width*0.5+coord;
275 Dist = (lindist>0) ? -lindist/Comp : 0;
276 signC= -1.0;
277 }
278 else
279 {
280 Dist = kInfinity;
281 }
282 // signC = sign<G4double>(Comp)
283 candidateNormal.exitNormal = ( signC * VecCartAxes[axis]);
284 candidateNormal.calculated = true;
285 candidateNormal.validConvex = true;
286 candidateNormal.exitSide =
287 (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
288 break;
289 case kPhi:
290 Dist = DistanceToOutPhi(localPoint,localDirection,width,candidateNormal);
291 // candidateNormal set in call
292 break;
293 case kRho:
294 Dist = DistanceToOutRad(localPoint,localDirection,width,offset,
295 replicaNo,candidateNormal);
296 // candidateNormal set in call
297 break;
298 default:
299 G4Exception("G4ReplicaNavigation::DistanceToOut()", "GeomNav0002",
300 FatalException, "Unknown axis!");
301 break;
302 }
303
304 arExitNormal= candidateNormal; // .exitNormal;
305
306 return Dist;
307 }
308
309 // ********************************************************************
310 // DistanceToOutPhi
311 // ********************************************************************
312 //
313 G4double
314 G4ReplicaNavigation::DistanceToOutPhi(const G4ThreeVector& localPoint,
315 const G4ThreeVector& localDirection,
316 const G4double width,
317 G4ExitNormal& foundNormal ) const
318 {
319 // Phi Intersection
320 // NOTE: width<=pi by definition
321 //
322 G4double sinSPhi = -2.0, cosSPhi = -2.0;
323 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
324 G4ExitNormal::ESide sidePhi = G4ExitNormal::kNull;
325 G4ThreeVector candidateNormal;
326
327 if ( (localPoint.x()!=0.0) || (localPoint.y()!=0.0) )
328 {
329 sinSPhi = std::sin(-width*0.5); // SIN of starting phi plane
330 cosSPhi = std::cos(width*0.5); // COS of starting phi plane
331
332 // pDist -ve when inside
333 //
334 pDistS = localPoint.x()*sinSPhi-localPoint.y()*cosSPhi;
335 // Start plane at phi= -S
336 pDistE = localPoint.x()*sinSPhi+localPoint.y()*cosSPhi;
337 // End plane at phi= +S
338
339 // Comp -ve when in direction of outwards normal
340 //
341 compS = -sinSPhi*localDirection.x()+cosSPhi*localDirection.y();
342 compE = -sinSPhi*localDirection.x()-cosSPhi*localDirection.y();
343
344 if ( (pDistS<=halfkCarTolerance)&&(pDistE<=halfkCarTolerance) )
345 {
346 // Inside both phi *full* planes
347 //
348 if ( compS<0 )
349 {
350 dist2 = pDistS/compS;
351 yi = localPoint.y()+dist2*localDirection.y();
352
353 // Check intersecting with correct half-plane (no -> no intersect)
354 //
355 if ( yi<=0 )
356 {
357 Dist = (pDistS<=-halfkCarTolerance) ? dist2 : 0;
358 sidePhi= G4ExitNormal::kSPhi; // tbc
359 }
360 else
361 {
362 Dist = kInfinity;
363 }
364 }
365 else
366 {
367 Dist = kInfinity;
368 }
369 if ( compE<0 )
370 {
371 dist2 = pDistE/compE;
372
373 // Only check further if < starting phi intersection
374 //
375 if ( dist2<Dist )
376 {
377 yi = localPoint.y()+dist2*localDirection.y();
378
379 // Check intersecting with correct half-plane
380 //
381 if ( yi>=0 )
382 {
383 // Leaving via ending phi
384 //
385 Dist = (pDistE<=-halfkCarTolerance) ? dist2 : 0;
386 sidePhi = G4ExitNormal::kEPhi;
387 }
388 }
389 }
390 }
391 else if ( (pDistS>halfkCarTolerance)&&(pDistE>halfkCarTolerance) )
392 {
393 // Outside both *full* phi planes
394 // if towards both >=0 then once inside will remain inside
395 //
396 Dist = ((compS>=0)&&(compE>=0)) ? kInfinity : 0;
397 }
398 else if ( (pDistS>halfkCarTolerance)&&(pDistE<=halfkCarTolerance) )
399 {
400 // Outside full starting plane, inside full ending plane
401 //
402 if ( compE<0 )
403 {
404 dist2 = pDistE/compE;
405 yi = localPoint.y()+dist2*localDirection.y();
406
407 // Check intersection in correct half-plane
408 // (if not -> remain in extent)
409 //
410 Dist = (yi>0) ? dist2 : kInfinity;
411 if( yi> 0 ) { sidePhi = G4ExitNormal::kEPhi; }
412 }
413 else // Leaving immediately by starting phi
414 {
415 Dist = kInfinity;
416 }
417 }
418 else
419 {
420 // Must be (pDistS<=halfkCarTolerance)&&(pDistE>halfkCarTolerance)
421 // Inside full starting plane, outside full ending plane
422 //
423 if ( compE>=0 )
424 {
425 if ( compS<0 )
426 {
427 dist2 = pDistS/compS;
428 yi = localPoint.y()+dist2*localDirection.y();
429
430 // Check intersection in correct half-plane
431 // (if not -> remain in extent)
432 //
433 Dist = (yi<0) ? dist2 : kInfinity;
434 if(yi<0) { sidePhi = G4ExitNormal::kSPhi; }
435 }
436 else
437 {
438 Dist = kInfinity;
439 }
440 }
441 else
442 {
443 // Leaving immediately by ending phi
444 //
445 Dist = 0.;
446 sidePhi = G4ExitNormal::kEPhi;
447 }
448 }
449 }
450 else
451 {
452 // On z axis + travel not || to z axis -> use direction vector
453 //
454 if( (std::fabs(localDirection.phi())<=width*0.5) )
455 {
456 Dist = kInfinity;
457 }
458 else
459 {
460 Dist = 0.;
461 sidePhi = G4ExitNormal::kMY;
462 }
463 }
464
465 if(sidePhi == G4ExitNormal::kSPhi )
466 {
467 candidateNormal = G4ThreeVector(sinSPhi,-cosSPhi,0.) ;
468 }
469 else if (sidePhi == G4ExitNormal::kEPhi)
470 {
471 candidateNormal = G4ThreeVector(sinSPhi,cosSPhi,0.) ;
472 }
473 else if (sidePhi == G4ExitNormal::kMY )
474 {
475 candidateNormal = G4ThreeVector(0., -1.0, 0.); // Split -S and +S 'phi'
476 }
477 foundNormal.calculated= (sidePhi != G4ExitNormal::kNull );
478 foundNormal.exitNormal= candidateNormal;
479
480 return Dist;
481 }
482
483 // ********************************************************************
484 // DistanceToOutRad
485 // ********************************************************************
486 //
487 G4double
488 G4ReplicaNavigation::DistanceToOutRad(const G4ThreeVector& localPoint,
489 const G4ThreeVector& localDirection,
490 const G4double width,
491 const G4double offset,
492 const G4int replicaNo,
493 G4ExitNormal& foundNormal ) const
494 {
495 G4double rmin, rmax, t1, t2, t3, deltaR;
496 G4double b, c, d2, srd;
497 G4ExitNormal::ESide sideR= G4ExitNormal::kNull;
498
499 //
500 // Radial Intersections
501 //
502
503 // Find intersction with cylinders at rmax/rmin
504 // Intersection point (xi,yi,zi) on line
505 // x=localPoint.x+t*localDirection.x etc.
506 //
507 // Intersects with x^2+y^2=R^2
508 //
509 // Hence (localDirection.x^2+localDirection.y^2)t^2+
510 // 2t(localPoint.x*localDirection.x+localPoint.y*localDirection.y)+
511 // localPoint.x^2+localPoint.y^2-R^2=0
512 //
513 // t1 t2 t3
514
515 rmin = replicaNo*width+offset;
516 rmax = (replicaNo+1)*width+offset;
517
518 t1 = 1.0-localDirection.z()*localDirection.z(); // since v normalised
519 t2 = localPoint.x()*localDirection.x()+localPoint.y()*localDirection.y();
520 t3 = localPoint.x()*localPoint.x()+localPoint.y()*localPoint.y();
521
522 if ( t1>0 ) // Check not parallel
523 {
524 // Calculate srd, r exit distance
525 //
526 if ( t2>=0 )
527 {
528 // Delta r not negative => leaving via rmax
529 //
530 deltaR = t3-rmax*rmax;
531
532 // NOTE: Should use
533 // rho-rmax<-halfkRadTolerance - [no sqrts for efficiency]
534 //
535 if ( deltaR<-halfkRadTolerance )
536 {
537 b = t2/t1;
538 c = deltaR/t1;
539 srd = -b+std::sqrt(b*b-c);
540 sideR = G4ExitNormal::kRMax;
541 }
542 else
543 {
544 // On tolerant boundary & heading outwards (or locally
545 // perpendicular to) outer radial surface -> leaving immediately
546 //
547 srd = 0;
548 sideR = G4ExitNormal::kRMax;
549 }
550 }
551 else
552 {
553 // Possible rmin intersection
554 //
555 if (rmin != 0.0)
556 {
557 deltaR = t3-rmin*rmin;
558 b = t2/t1;
559 c = deltaR/t1;
560 d2 = b*b-c;
561 if ( d2>=0 )
562 {
563 // Leaving via rmin
564 // NOTE: Should use
565 // rho-rmin>halfkRadTolerance - [no sqrts for efficiency]
566 //
567 srd = (deltaR>halfkRadTolerance) ? -b-std::sqrt(d2) : 0.0;
568 // Is the following more accurate ?
569 // srd = (deltaR>halfkRadTolerance) ? c/( -b - std::sqrt(d2)) : 0.0;
570 sideR = G4ExitNormal::kRMin;
571 }
572 else
573 {
574 // No rmin intersect -> must be rmax intersect
575 //
576 deltaR = t3-rmax*rmax;
577 c = deltaR/t1;
578 d2 = b*b-c;
579 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
580 sideR = G4ExitNormal::kRMax;
581 }
582 }
583 else
584 {
585 // No rmin intersect -> must be rmax intersect
586 //
587 deltaR = t3-rmax*rmax;
588 b = t2/t1;
589 c = deltaR/t1;
590 d2 = b*b-c;
591 srd = (d2 < 0.) ? 0.0 : -b+std::sqrt(d2);
592 sideR= G4ExitNormal::kRMax;
593 }
594 }
595 }
596 else
597 {
598 srd =kInfinity;
599 sideR = G4ExitNormal::kNull;
600 }
601
602 if( sideR != G4ExitNormal::kNull ) // if ((side == kRMax) || (side==kRMin))
603 {
604 // Note: returned vector not explicitly normalised
605 // (divided by fRMax for unit vector)
606
607 G4double xi, yi;
608 xi = localPoint.x() + srd*localDirection.x();
609 yi = localPoint.y() + srd*localDirection.y();
610 G4ThreeVector normalR = G4ThreeVector(xi,yi,0.0);
611
612 if( sideR == G4ExitNormal::kRMax )
613 {
614 normalR *= 1.0/rmax;
615 }
616 else
617 {
618 normalR *= (-1.0)/rmin;
619 }
620 foundNormal.exitNormal= normalR;
621 foundNormal.calculated= true;
622 foundNormal.validConvex = (sideR == G4ExitNormal::kRMax);
623 foundNormal.exitSide = sideR;
624 }
625 else
626 {
627 foundNormal.calculated = false;
628 }
629
630 return srd;
631 }
632
633 // ********************************************************************
634 // ComputeTransformation
635 //
636 // Setup transformation and transform point into local system
637 // ********************************************************************
638 //
639 void
640 G4ReplicaNavigation::ComputeTransformation(const G4int replicaNo,
641 G4VPhysicalVolume* pVol,
642 G4ThreeVector& point) const
643 {
644 G4double val,cosv,sinv,tmpx,tmpy;
645
646 // Replication data
647 //
648 EAxis axis;
649 G4int nReplicas;
650 G4double width,offset;
651 G4bool consuming;
652
653 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
654
655 switch (axis)
656 {
657 case kXAxis:
658 val = -width*0.5*(nReplicas-1)+width*replicaNo;
659 pVol->SetTranslation(G4ThreeVector(val,0,0));
660 point.setX(point.x()-val);
661 break;
662 case kYAxis:
663 val = -width*0.5*(nReplicas-1)+width*replicaNo;
664 pVol->SetTranslation(G4ThreeVector(0,val,0));
665 point.setY(point.y()-val);
666 break;
667 case kZAxis:
668 val = -width*0.5*(nReplicas-1)+width*replicaNo;
669 pVol->SetTranslation(G4ThreeVector(0,0,val));
670 point.setZ(point.z()-val);
671 break;
672 case kPhi:
673 val = -(offset+width*(replicaNo+0.5));
674 SetPhiTransformation(val,pVol);
675 cosv = std::cos(val);
676 sinv = std::sin(val);
677 tmpx = point.x()*cosv-point.y()*sinv;
678 tmpy = point.x()*sinv+point.y()*cosv;
679 point.setY(tmpy);
680 point.setX(tmpx);
681 break;
682 case kRho:
683 // No setup required for radial case
684 default:
685 break;
686 }
687 }
688
689 // ********************************************************************
690 // ComputeTransformation
691 //
692 // Setup transformation into local system
693 // ********************************************************************
694 //
695 void
696 G4ReplicaNavigation::ComputeTransformation(const G4int replicaNo,
697 G4VPhysicalVolume* pVol) const
698 {
699 G4double val;
700
701 // Replication data
702 //
703 EAxis axis;
704 G4int nReplicas;
705 G4double width, offset;
706 G4bool consuming;
707
708 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
709
710 switch (axis)
711 {
712 case kXAxis:
713 val = -width*0.5*(nReplicas-1)+width*replicaNo;
714 pVol->SetTranslation(G4ThreeVector(val,0,0));
715 break;
716 case kYAxis:
717 val = -width*0.5*(nReplicas-1)+width*replicaNo;
718 pVol->SetTranslation(G4ThreeVector(0,val,0));
719 break;
720 case kZAxis:
721 val = -width*0.5*(nReplicas-1)+width*replicaNo;
722 pVol->SetTranslation(G4ThreeVector(0,0,val));
723 break;
724 case kPhi:
725 val = -(offset+width*(replicaNo+0.5));
726 SetPhiTransformation(val,pVol);
727 break;
728 case kRho:
729 // No setup required for radial case
730 default:
731 break;
732 }
733 }
734
735 // ********************************************************************
736 // ComputeStep
737 // ********************************************************************
738 //
739 G4double
740 G4ReplicaNavigation::ComputeStep(const G4ThreeVector& globalPoint,
741 const G4ThreeVector& globalDirection,
742 const G4ThreeVector& localPoint,
743 const G4ThreeVector& localDirection,
744 const G4double currentProposedStepLength,
745 G4double& newSafety,
746 G4NavigationHistory &history,
747 // std::pair<G4bool,G4bool>& validAndCalculated
748 G4bool& validExitNormal,
749 G4bool& calculatedExitNormal,
750 G4ThreeVector& exitNormalVector,
751 G4bool& exiting,
752 G4bool& entering,
753 G4VPhysicalVolume* (*pBlockedPhysical),
754 G4int& blockedReplicaNo )
755 {
756 G4VPhysicalVolume *repPhysical, *motherPhysical;
757 G4VPhysicalVolume *samplePhysical, *blockedExitedVol = nullptr;
758 G4LogicalVolume *repLogical;
759 G4VSolid *motherSolid;
760 G4ThreeVector repPoint, repDirection, sampleDirection;
761 G4double ourStep=currentProposedStepLength;
762 G4double ourSafety=kInfinity;
763 G4double sampleStep, sampleSafety, motherStep, motherSafety;
764 G4long localNoDaughters, sampleNo;
765 G4int depth;
766 G4ExitNormal exitNormalStc;
767 // G4int depthDeterminingStep= -1; // Useful only for debugging - for now
768
769 calculatedExitNormal= false;
770
771 // Exiting normal optimisation
772 //
773 if ( exiting&&validExitNormal )
774 {
775 if ( localDirection.dot(exitNormalVector)>=kMinExitingNormalCosine )
776 {
777 // Block exited daughter volume
778 //
779 blockedExitedVol = *pBlockedPhysical;
780 ourSafety = 0;
781 }
782 }
783 exiting = false;
784 entering = false;
785
786 repPhysical = history.GetTopVolume();
787 repLogical = repPhysical->GetLogicalVolume();
788
789 //
790 // Compute intersection with replica boundaries & replica safety
791 //
792
793 sampleSafety = DistanceToOut(repPhysical,
794 history.GetTopReplicaNo(),
795 localPoint);
796 G4ExitNormal normalOutStc;
797 const auto topDepth= (G4int)history.GetDepth();
798
799 ourSafety = std::min( ourSafety, sampleSafety);
800
801 if ( sampleSafety<ourStep )
802 {
803
804 sampleStep = DistanceToOut(repPhysical,
805 history.GetTopReplicaNo(),
806 localPoint,
807 localDirection,
808 normalOutStc);
809 if ( sampleStep<ourStep )
810 {
811 ourStep = sampleStep;
812 exiting = true;
813 validExitNormal = normalOutStc.validConvex; // false; -> Old,Conservative
814
815 exitNormalStc = normalOutStc;
816 exitNormalStc.exitNormal =
817 history.GetTopTransform().InverseTransformAxis(normalOutStc.exitNormal);
818 calculatedExitNormal = true;
819 }
820 }
821 const G4int secondDepth = topDepth;
822 depth = secondDepth;
823
824 // Loop checking, 07.10.2016, JA -- Need to add: assert(depth>0)
825 while ( history.GetVolumeType(depth)==kReplica )
826 {
827 const G4AffineTransform& GlobalToLocal = history.GetTransform(depth);
828 repPoint = GlobalToLocal.TransformPoint(globalPoint);
829 // repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
830
831 sampleSafety = DistanceToOut(history.GetVolume(depth),
832 history.GetReplicaNo(depth),
833 repPoint);
834 if ( sampleSafety < ourSafety )
835 {
836 ourSafety = sampleSafety;
837 }
838 if ( sampleSafety < ourStep )
839 {
840 G4ThreeVector newLocalDirection =
841 GlobalToLocal.TransformAxis(globalDirection);
842 sampleStep = DistanceToOut(history.GetVolume(depth),
843 history.GetReplicaNo(depth),
844 repPoint,
845 newLocalDirection,
846 normalOutStc);
847 if ( sampleStep < ourStep )
848 {
849 ourStep = sampleStep;
850 exiting = true;
851
852 // As step is limited by this level, must set Exit Normal
853 //
854 G4ThreeVector localExitNorm = normalOutStc.exitNormal;
855 G4ThreeVector globalExitNorm =
856 GlobalToLocal.InverseTransformAxis(localExitNorm);
857
858 exitNormalStc = normalOutStc; // Normal, convex, calculated, side
859 exitNormalStc.exitNormal = globalExitNorm;
860 calculatedExitNormal = true;
861 }
862 }
863 depth--;
864 }
865
866 // Compute mother safety & intersection
867 //
868 G4ThreeVector exitVectorMother;
869 G4bool exitConvex = false; // Value obtained in DistanceToOut(p,v) call
870 G4ExitNormal motherNormalStc;
871
872 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
873 motherPhysical = history.GetVolume(depth);
874 motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
875 motherSafety = motherSolid->DistanceToOut(repPoint);
876 repDirection = history.GetTransform(depth).TransformAxis(globalDirection);
877
878 motherStep = motherSolid->DistanceToOut(repPoint,repDirection,true,
879 &exitConvex,&exitVectorMother);
880 if( exitConvex )
881 {
882 motherNormalStc = G4ExitNormal( exitVectorMother, true, false,
883 G4ExitNormal::kMother);
884 calculatedExitNormal = true;
885 }
886 const G4AffineTransform& globalToLocalTop = history.GetTopTransform();
887
888 G4bool motherDeterminedStep = (motherStep<ourStep);
889
890 if( (!exitConvex) && motherDeterminedStep )
891 {
892 exitVectorMother = motherSolid->SurfaceNormal( repPoint );
893 motherNormalStc = G4ExitNormal( exitVectorMother, true, false,
894 G4ExitNormal::kMother);
895 // CalculatedExitNormal -> true;
896 // Convex -> false: do not know value
897 // ExitSide -> kMother (or kNull)
898
899 calculatedExitNormal = true;
900 }
901 if( motherDeterminedStep )
902 {
903 G4ThreeVector globalExitNormalTop =
904 globalToLocalTop.InverseTransformAxis(exitVectorMother);
905
906 exitNormalStc = motherNormalStc;
907 exitNormalStc.exitNormal = globalExitNormalTop;
908 }
909
910 // Push in principle no longer necessary. G4Navigator now takes care of ...
911 // Removing this however may cause additional almost-zero steps and generate
912 // warnings for pushed particles from G4Navigator, particularly for the case
913 // of 3D replicas (Cartesian or combined Radial/Phi cases).
914 // Requires further investigation and eventually reimplementation of
915 // LevelLocate() to take into account point and direction ...
916 //
917 if( ourStep<fMinStep )
918 {
919 ourStep = 2*kCarTolerance;
920 }
921
922 if ( motherSafety<ourSafety )
923 {
924 ourSafety = motherSafety;
925 }
926
927 #ifdef G4VERBOSE
928 if ( fCheck )
929 {
930 if( motherSolid->Inside(localPoint)==kOutside )
931 {
932 std::ostringstream message;
933 message << "Point outside volume !" << G4endl
934 << " Point " << localPoint
935 << " is outside current volume " << motherPhysical->GetName()
936 << G4endl;
937 G4double estDistToSolid= motherSolid->DistanceToIn(localPoint);
938 message << " Estimated isotropic distance to solid (distToIn)= "
939 << estDistToSolid << G4endl;
940 if( estDistToSolid > 100.0 * kCarTolerance )
941 {
942 motherSolid->DumpInfo();
943 G4Exception("G4ReplicaNavigation::ComputeStep()",
944 "GeomNav0003", FatalException, message,
945 "Point is far outside Current Volume !" );
946 }
947 else
948 {
949 G4Exception("G4ReplicaNavigation::ComputeStep()",
950 "GeomNav1002", JustWarning, message,
951 "Point is a little outside Current Volume.");
952 }
953 }
954 }
955 #endif
956
957 // Comparison of steps may need precision protection
958 //
959 #if 1
960 if( motherDeterminedStep )
961 {
962 ourStep = motherStep;
963 exiting = true;
964 }
965
966 // Transform it to the Grand-Mother Reference Frame (current convention)
967 //
968 if ( calculatedExitNormal )
969 {
970 if ( motherDeterminedStep )
971 {
972 exitNormalVector = motherNormalStc.exitNormal;
973 }
974 else
975 {
976 G4ThreeVector exitNormalGlobal = exitNormalStc.exitNormal;
977 exitNormalVector = globalToLocalTop.TransformAxis(exitNormalGlobal);
978 // exitNormalVector= globalToLocal2nd.TransformAxis(exitNormalGlobal);
979 // Alt Make it in one go to Grand-Mother, avoiding transform below
980 }
981 // Transform to Grand-mother reference frame
982 const G4RotationMatrix* rot = motherPhysical->GetRotation();
983 if ( rot != nullptr )
984 {
985 exitNormalVector *= rot->inverse();
986 }
987
988 }
989 else
990 {
991 validExitNormal = false;
992 }
993
994 #else
995 if ( motherSafety<=ourStep )
996 {
997 if ( motherStep<=ourStep )
998 {
999 ourStep = motherStep;
1000 exiting = true;
1001 if ( validExitNormal )
1002 {
1003 const G4RotationMatrix* rot = motherPhysical->GetRotation();
1004 if ( rot )
1005 {
1006 exitNormal *= rot->inverse();
1007 }
1008 }
1009 }
1010 else
1011 {
1012 validExitNormal = false;
1013 // calculatedExitNormal= false;
1014 }
1015 }
1016 #endif
1017
1018
1019 G4bool daughterDeterminedStep = false;
1020 G4ThreeVector daughtNormRepCrd;
1021 // Exit normal of daughter transformed to
1022 // the coordinate system of Replica (i.e. last depth)
1023
1024 //
1025 // Compute daughter safeties & intersections
1026 //
1027 localNoDaughters = repLogical->GetNoDaughters();
1028 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1029 {
1030 samplePhysical = repLogical->GetDaughter((G4int)sampleNo);
1031 if ( samplePhysical!=blockedExitedVol )
1032 {
1033 G4ThreeVector localExitNorm;
1034 G4ThreeVector normReplicaCoord;
1035
1036 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
1037 samplePhysical->GetTranslation());
1038 sampleTf.Invert();
1039 const G4ThreeVector samplePoint =
1040 sampleTf.TransformPoint(localPoint);
1041 const G4VSolid* sampleSolid =
1042 samplePhysical->GetLogicalVolume()->GetSolid();
1043 const G4double sampleSafetyDistance =
1044 sampleSolid->DistanceToIn(samplePoint);
1045 if ( sampleSafetyDistance<ourSafety )
1046 {
1047 ourSafety = sampleSafetyDistance;
1048 }
1049 if ( sampleSafetyDistance<=ourStep )
1050 {
1051 sampleDirection = sampleTf.TransformAxis(localDirection);
1052 const G4double sampleStepDistance =
1053 sampleSolid->DistanceToIn(samplePoint,sampleDirection);
1054 if ( sampleStepDistance<=ourStep )
1055 {
1056 daughterDeterminedStep = true;
1057
1058 ourStep = sampleStepDistance;
1059 entering = true;
1060 exiting = false;
1061 *pBlockedPhysical = samplePhysical;
1062 blockedReplicaNo = (G4int)sampleNo;
1063
1064 #ifdef DAUGHTER_NORMAL_ALSO
1065 // This norm can be calculated later, if needed daughter is available
1066 localExitNorm = sampleSolid->SurfaceNormal(samplePoint);
1067 daughtNormRepCrd = sampleTf.InverseTransformAxis(localExitNorm);
1068 #endif
1069
1070 #ifdef G4VERBOSE
1071 // Check to see that the resulting point is indeed in/on volume.
1072 // This check could eventually be made only for successful candidate.
1073
1074 if ( ( fCheck ) && ( sampleStepDistance < kInfinity ) )
1075 {
1076 G4ThreeVector intersectionPoint;
1077 intersectionPoint = samplePoint
1078 + sampleStepDistance * sampleDirection;
1079 EInside insideIntPt = sampleSolid->Inside(intersectionPoint);
1080 if ( insideIntPt != kSurface )
1081 {
1082 G4long oldcoutPrec = G4cout.precision(16);
1083 std::ostringstream message;
1084 message << "Navigator gets conflicting response from Solid."
1085 << G4endl
1086 << " Inaccurate DistanceToIn for solid "
1087 << sampleSolid->GetName() << G4endl
1088 << " Solid gave DistanceToIn = "
1089 << sampleStepDistance << " yet returns " ;
1090 if ( insideIntPt == kInside ) {
1091 message << "-kInside-";
1092 } else if ( insideIntPt == kOutside ) {
1093 message << "-kOutside-";
1094 } else {
1095 message << "-kSurface-";
1096 }
1097 message << " for this point !" << G4endl
1098 << " Point = " << intersectionPoint << G4endl;
1099 if ( insideIntPt != kInside ) {
1100 message << " DistanceToIn(p) = "
1101 << sampleSolid->DistanceToIn(intersectionPoint)
1102 << G4endl;
1103 }
1104 if ( insideIntPt != kOutside ) {
1105 message << " DistanceToOut(p) = "
1106 << sampleSolid->DistanceToOut(intersectionPoint);
1107 }
1108 G4Exception("G4ReplicaNavigation::ComputeStep()",
1109 "GeomNav1002", JustWarning, message);
1110 G4cout.precision(oldcoutPrec);
1111 }
1112 }
1113 #endif
1114 }
1115 }
1116 }
1117 }
1118
1119 calculatedExitNormal &= (!daughterDeterminedStep);
1120
1121 #ifdef DAUGHTER_NORMAL_ALSO
1122 if( daughterDeterminedStep )
1123 {
1124 // G4ThreeVector daughtNormGlobal =
1125 // GlobalToLastDepth.Inverse().TransformAxis(daughtNormRepCrd);
1126 // ==> Can calculate it, but have no way to transmit it to caller (for now)
1127
1128 exitNormalVector = globalToLocalTop.InverseTransformAxis(daughtNormGlobal);
1129 validExitNormal = false; // Entering daughter - never convex for parent
1130
1131 calculatedExitNormal = true;
1132 }
1133 // calculatedExitNormal= true; // Force it to true -- dubious
1134 #endif
1135
1136 newSafety = ourSafety;
1137 return ourStep;
1138 }
1139
1140 // ********************************************************************
1141 // ComputeSafety
1142 //
1143 // Compute the isotropic distance to current volume's boundaries
1144 // and to daughter volumes.
1145 // ********************************************************************
1146 //
1147 G4double
1148 G4ReplicaNavigation::ComputeSafety(const G4ThreeVector& globalPoint,
1149 const G4ThreeVector& localPoint,
1150 const G4NavigationHistory& history,
1151 const G4double ) const
1152 {
1153 G4VPhysicalVolume *repPhysical, *motherPhysical;
1154 G4VPhysicalVolume *samplePhysical, *blockedExitedVol = nullptr;
1155 G4LogicalVolume *repLogical;
1156 G4VSolid *motherSolid;
1157 G4ThreeVector repPoint;
1158 G4double ourSafety = kInfinity;
1159 G4double sampleSafety;
1160 G4long localNoDaughters, sampleNo;
1161 G4int depth;
1162
1163 repPhysical = history.GetTopVolume();
1164 repLogical = repPhysical->GetLogicalVolume();
1165
1166 //
1167 // Compute intersection with replica boundaries & replica safety
1168 //
1169
1170 sampleSafety = DistanceToOut(history.GetTopVolume(),
1171 history.GetTopReplicaNo(),
1172 localPoint);
1173 if ( sampleSafety<ourSafety )
1174 {
1175 ourSafety = sampleSafety;
1176 }
1177
1178 depth = (G4int)history.GetDepth()-1;
1179
1180 // Loop checking, 07.10.2016, JA -- need to add: assert(depth>0)
1181 while ( history.GetVolumeType(depth)==kReplica )
1182 {
1183 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1184 sampleSafety = DistanceToOut(history.GetVolume(depth),
1185 history.GetReplicaNo(depth),
1186 repPoint);
1187 if ( sampleSafety<ourSafety )
1188 {
1189 ourSafety = sampleSafety;
1190 }
1191 depth--;
1192 }
1193
1194 // Compute mother safety & intersection
1195 //
1196 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1197 motherPhysical = history.GetVolume(depth);
1198 motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
1199 sampleSafety = motherSolid->DistanceToOut(repPoint);
1200
1201 if ( sampleSafety<ourSafety )
1202 {
1203 ourSafety = sampleSafety;
1204 }
1205
1206 // Compute daughter safeties & intersections
1207 //
1208 localNoDaughters = repLogical->GetNoDaughters();
1209 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1210 {
1211 samplePhysical = repLogical->GetDaughter((G4int)sampleNo);
1212 if ( samplePhysical!=blockedExitedVol )
1213 {
1214 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
1215 samplePhysical->GetTranslation());
1216 sampleTf.Invert();
1217 const G4ThreeVector samplePoint =
1218 sampleTf.TransformPoint(localPoint);
1219 const G4VSolid *sampleSolid =
1220 samplePhysical->GetLogicalVolume()->GetSolid();
1221 const G4double sampleSafetyDistance =
1222 sampleSolid->DistanceToIn(samplePoint);
1223 if ( sampleSafetyDistance<ourSafety )
1224 {
1225 ourSafety = sampleSafetyDistance;
1226 }
1227 }
1228 }
1229 return ourSafety;
1230 }
1231
1232 // ********************************************************************
1233 // BackLocate
1234 // ********************************************************************
1235 //
1236 EInside
1237 G4ReplicaNavigation::BackLocate(G4NavigationHistory& history,
1238 const G4ThreeVector& globalPoint,
1239 G4ThreeVector& localPoint,
1240 const G4bool& exiting,
1241 G4bool& notKnownInside ) const
1242 {
1243 G4VPhysicalVolume *pNRMother = nullptr;
1244 G4VSolid *motherSolid;
1245 G4ThreeVector repPoint, goodPoint;
1246 G4int mdepth, depth, cdepth;
1247 EInside insideCode;
1248
1249 cdepth = (G4int)history.GetDepth();
1250
1251 // Find non replicated mother
1252 //
1253 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1254 {
1255 if ( history.GetVolumeType(mdepth)!=kReplica )
1256 {
1257 pNRMother = history.GetVolume(mdepth);
1258 break;
1259 }
1260 }
1261
1262 if( pNRMother == nullptr )
1263 {
1264 // All the tree of mother volumes were Replicas.
1265 // This is an error, as the World volume must be a Placement
1266 //
1267 G4Exception("G4ReplicaNavigation::BackLocate()", "GeomNav0002",
1268 FatalException, "The World volume must be a Placement!");
1269 return kInside;
1270 }
1271
1272 motherSolid = pNRMother->GetLogicalVolume()->GetSolid();
1273 goodPoint = history.GetTransform(mdepth).TransformPoint(globalPoint);
1274 insideCode = motherSolid->Inside(goodPoint);
1275 if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1276 {
1277 // Outside mother -> back up to mother level
1278 // Locate.. in Navigator will back up one more level
1279 // localPoint not required
1280 //
1281 history.BackLevel(cdepth-mdepth);
1282 // localPoint = goodPoint;
1283 }
1284 else
1285 {
1286 notKnownInside = false;
1287
1288 // Still within replications
1289 // Check down: if on outside stop at this level
1290 //
1291 for ( depth=mdepth+1; depth<cdepth; ++depth)
1292 {
1293 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1294 insideCode = Inside(history.GetVolume(depth),
1295 history.GetReplicaNo(depth),
1296 repPoint);
1297 if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1298 {
1299 localPoint = goodPoint;
1300 history.BackLevel(cdepth-depth);
1301 return insideCode;
1302 }
1303 goodPoint = repPoint;
1304 }
1305 localPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1306 insideCode = Inside(history.GetVolume(depth),
1307 history.GetReplicaNo(depth),
1308 localPoint);
1309 // If outside level, set localPoint = coordinates in reference system
1310 // of *previous* level - location code in navigator will back up one
1311 // level [And also manage blocking]
1312 //
1313 if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1314 {
1315 localPoint = goodPoint;
1316 }
1317 }
1318 return insideCode;
1319 }
1320