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