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