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