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