Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 100372 2016-10-19 10:37:52Z 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  halfkCarTolerance = kCarTolerance*0.5;
64  halfkRadTolerance = kRadTolerance*0.5;
65  halfkAngTolerance = kAngTolerance*0.5;
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
328 G4ReplicaNavigation::DistanceToOutPhi(const G4ThreeVector &localPoint,
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<=halfkCarTolerance)&&(pDistE<=halfkCarTolerance) )
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>halfkCarTolerance)&&(pDistE>halfkCarTolerance) )
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>halfkCarTolerance)&&(pDistE<=halfkCarTolerance) )
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<=halfkCarTolerance)&&(pDistE>halfkCarTolerance)
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
502 G4ReplicaNavigation::DistanceToOutRad(const G4ThreeVector &localPoint,
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 
621  G4double xi, yi;
622  xi = localPoint.x() + srd*localDirection.x();
623  yi = localPoint.y() + srd*localDirection.y();
624  G4ThreeVector normalR = G4ThreeVector(xi,yi,0.0);
625 
626  if( sideR == G4ExitNormal::kRMax )
627  {
628  normalR *= 1.0/rmax;
629  }
630  else
631  {
632  normalR *= (-1.0)/rmin;
633  }
634  foundNormal.exitNormal= normalR;
635  foundNormal.calculated= true;
636  foundNormal.validConvex = (sideR == G4ExitNormal::kRMax);
637  foundNormal.exitSide = sideR;
638  }
639  else
640  {
641  foundNormal.calculated= false;
642  }
643 
644  return srd;
645 }
646 
647 // ********************************************************************
648 // ComputeTransformation
649 //
650 // Setup transformation and transform point into local system
651 // ********************************************************************
652 //
653 void
655  G4VPhysicalVolume* pVol,
656  G4ThreeVector& point) const
657 {
658  G4double val,cosv,sinv,tmpx,tmpy;
659 
660  // Replication data
661  //
662  EAxis axis;
663  G4int nReplicas;
664  G4double width,offset;
665  G4bool consuming;
666 
667  pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
668 
669  switch (axis)
670  {
671  case kXAxis:
672  val = -width*0.5*(nReplicas-1)+width*replicaNo;
673  pVol->SetTranslation(G4ThreeVector(val,0,0));
674  point.setX(point.x()-val);
675  break;
676  case kYAxis:
677  val = -width*0.5*(nReplicas-1)+width*replicaNo;
678  pVol->SetTranslation(G4ThreeVector(0,val,0));
679  point.setY(point.y()-val);
680  break;
681  case kZAxis:
682  val = -width*0.5*(nReplicas-1)+width*replicaNo;
683  pVol->SetTranslation(G4ThreeVector(0,0,val));
684  point.setZ(point.z()-val);
685  break;
686  case kPhi:
687  val = -(offset+width*(replicaNo+0.5));
688  SetPhiTransformation(val,pVol);
689  cosv = std::cos(val);
690  sinv = std::sin(val);
691  tmpx = point.x()*cosv-point.y()*sinv;
692  tmpy = point.x()*sinv+point.y()*cosv;
693  point.setY(tmpy);
694  point.setX(tmpx);
695  break;
696  case kRho:
697  // No setup required for radial case
698  default:
699  break;
700  }
701 }
702 
703 // ********************************************************************
704 // ComputeTransformation
705 //
706 // Setup transformation into local system
707 // ********************************************************************
708 //
709 void
711  G4VPhysicalVolume* pVol) const
712 {
713  G4double val;
714 
715  // Replication data
716  //
717  EAxis axis;
718  G4int nReplicas;
719  G4double width, offset;
720  G4bool consuming;
721 
722  pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
723 
724  switch (axis)
725  {
726  case kXAxis:
727  val = -width*0.5*(nReplicas-1)+width*replicaNo;
728  pVol->SetTranslation(G4ThreeVector(val,0,0));
729  break;
730  case kYAxis:
731  val = -width*0.5*(nReplicas-1)+width*replicaNo;
732  pVol->SetTranslation(G4ThreeVector(0,val,0));
733  break;
734  case kZAxis:
735  val = -width*0.5*(nReplicas-1)+width*replicaNo;
736  pVol->SetTranslation(G4ThreeVector(0,0,val));
737  break;
738  case kPhi:
739  val = -(offset+width*(replicaNo+0.5));
740  SetPhiTransformation(val,pVol);
741  break;
742  case kRho:
743  // No setup required for radial case
744  default:
745  break;
746  }
747 }
748 
749 // ********************************************************************
750 // ComputeStep
751 // ********************************************************************
752 //
753 G4double
755  const G4ThreeVector &globalDirection,
756  const G4ThreeVector &localPoint,
757  const G4ThreeVector &localDirection,
758  const G4double currentProposedStepLength,
759  G4double &newSafety,
760  G4NavigationHistory &history,
761  // std::pair<G4bool,G4bool> &validAndCalculated
762  G4bool &validExitNormal,
763  G4bool &calculatedExitNormal,
764  G4ThreeVector &exitNormalVector,
765  G4bool &exiting,
766  G4bool &entering,
767  G4VPhysicalVolume *(*pBlockedPhysical),
768  G4int &blockedReplicaNo )
769 {
770  G4VPhysicalVolume *repPhysical, *motherPhysical;
771  G4VPhysicalVolume *samplePhysical, *blockedExitedVol=0;
772  G4LogicalVolume *repLogical;
773  G4VSolid *motherSolid;
774  G4ThreeVector repPoint, repDirection, sampleDirection;
775  G4double ourStep=currentProposedStepLength;
776  G4double ourSafety=kInfinity;
777  G4double sampleStep, sampleSafety, motherStep, motherSafety;
778  G4int localNoDaughters, sampleNo;
779  G4int depth;
780  G4ExitNormal exitNormalStc;
781  // G4int depthDeterminingStep= -1; // Useful only for debugging - for now
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  // Loop checking, 07.10.2016, J.Apostolakis -- Need to add: assert(depth>0)
839  while ( history.GetVolumeType(depth)==kReplica )
840  {
841  const G4AffineTransform& GlobalToLocal= history.GetTransform(depth);
842  repPoint = GlobalToLocal.TransformPoint(globalPoint);
843  // repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
844 
845  sampleSafety = DistanceToOut(history.GetVolume(depth),
846  history.GetReplicaNo(depth),
847  repPoint);
848  if ( sampleSafety < ourSafety )
849  {
850  ourSafety = sampleSafety;
851  }
852  if ( sampleSafety < ourStep )
853  {
854  G4ThreeVector newLocalDirection =
855  GlobalToLocal.TransformAxis(globalDirection);
856  sampleStep = DistanceToOut(history.GetVolume(depth),
857  history.GetReplicaNo(depth),
858  repPoint,
859  newLocalDirection,
860  normalOutStc);
861  if ( sampleStep < ourStep )
862  {
863  ourStep = sampleStep;
864  exiting = true;
865 
866  // As step is limited by this level, must set Exit Normal
867  //
868  G4ThreeVector localExitNorm= normalOutStc.exitNormal;
869  G4ThreeVector globalExitNorm=
870  GlobalToLocal.Inverse().TransformAxis(localExitNorm);
871 
872  exitNormalStc= normalOutStc; // Normal, convex, calculated, side
873  exitNormalStc.exitNormal= globalExitNorm;
874  calculatedExitNormal= true;
875  }
876  }
877  depth--;
878  }
879 
880  // Compute mother safety & intersection
881  //
882  G4ThreeVector exitVectorMother;
883  G4bool exitConvex= false; // Value obtained in DistanceToOut(p,v) call
884  G4ExitNormal motherNormalStc;
885 
886  repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
887  motherPhysical = history.GetVolume(depth);
888  motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
889  motherSafety = motherSolid->DistanceToOut(repPoint);
890  repDirection = history.GetTransform(depth).TransformAxis(globalDirection);
891 
892  motherStep = motherSolid->DistanceToOut(repPoint,repDirection,true,
893  &exitConvex,&exitVectorMother);
894  if( exitConvex )
895  {
896  motherNormalStc = G4ExitNormal( exitVectorMother, true, false,
898  calculatedExitNormal= true;
899  }
900  const G4AffineTransform& globalToLocalTop = history.GetTopTransform();
901 
902  G4bool motherDeterminedStep= (motherStep<ourStep);
903 
904  if( (!exitConvex) && motherDeterminedStep )
905  {
906  exitVectorMother= motherSolid->SurfaceNormal( repPoint );
907  motherNormalStc= G4ExitNormal( exitVectorMother, true, false,
909  // CalculatedExitNormal -> true;
910  // Convex -> false: do not know value
911  // ExitSide -> kMother (or kNull)
912 
913  calculatedExitNormal= true;
914  }
915  if( motherDeterminedStep)
916  {
917  G4ThreeVector globalExitNormalTop=
918  globalToLocalTop.Inverse().TransformAxis(exitVectorMother);
919 
920  exitNormalStc= motherNormalStc;
921  exitNormalStc.exitNormal= globalExitNormalTop;
922  }
923 
924  // Push in principle no longer necessary. G4Navigator now takes care of ...
925  // Removing this however may cause additional almost-zero steps and generate
926  // warnings for pushed particles from G4Navigator, particularly for the case
927  // of 3D replicas (Cartesian or combined Radial/Phi cases).
928  // Requires further investigation and eventually reimplementation of
929  // LevelLocate() to take into account point and direction ...
930  //
931  if ( ( (ourStep<fMinStep) && (sampleSafety<halfkCarTolerance) )
932  && ( repLogical->GetSolid()->Inside(localPoint)==kSurface ) )
933  {
934  ourStep = 100*kCarTolerance;
935  }
936 
937  if ( motherSafety<ourSafety )
938  {
939  ourSafety = motherSafety;
940  }
941 
942 #ifdef G4VERBOSE
943  if ( fCheck )
944  {
945  if( motherSolid->Inside(localPoint)==kOutside )
946  {
947  std::ostringstream message;
948  message << "Point outside volume !" << G4endl
949  << " Point " << localPoint
950  << " is outside current volume " << motherPhysical->GetName()
951  << G4endl;
952  G4double estDistToSolid= motherSolid->DistanceToIn(localPoint);
953  message << " Estimated isotropic distance to solid (distToIn)= "
954  << estDistToSolid << G4endl;
955  if( estDistToSolid > 100.0 * kCarTolerance )
956  {
957  motherSolid->DumpInfo();
958  G4Exception("G4ReplicaNavigation::ComputeStep()",
959  "GeomNav0003", FatalException, message,
960  "Point is far outside Current Volume !" );
961  }
962  else
963  G4Exception("G4ReplicaNavigation::ComputeStep()",
964  "GeomNav1002", JustWarning, message,
965  "Point is a little outside Current Volume.");
966  }
967  }
968 #endif
969 
970  // Comparison of steps may need precision protection
971  //
972 #if 1
973  if( motherDeterminedStep)
974  {
975  ourStep = motherStep;
976  exiting = true;
977  }
978 
979  // Transform it to the Grand-Mother Reference Frame (current convention)
980  //
981  if ( calculatedExitNormal )
982  {
983  if ( motherDeterminedStep )
984  {
985  exitNormalVector= motherNormalStc.exitNormal;
986  }
987  else
988  {
989  G4ThreeVector exitNormalGlobal= exitNormalStc.exitNormal;
990  exitNormalVector= globalToLocalTop.TransformAxis(exitNormalGlobal);
991  // exitNormalVector= globalToLocal2nd.TransformAxis(exitNormalGlobal);
992  // Alt Make it in one go to Grand-Mother, avoiding transform below
993  }
994  // Transform to Grand-mother reference frame
995  const G4RotationMatrix* rot = motherPhysical->GetRotation();
996  if ( rot )
997  {
998  exitNormalVector *= rot->inverse();
999  }
1000 
1001  }
1002  else
1003  {
1004  validExitNormal = false;
1005  }
1006 
1007 #else
1008  if ( motherSafety<=ourStep )
1009  {
1010  if ( motherStep<=ourStep )
1011  {
1012  ourStep = motherStep;
1013  exiting = true;
1014  if ( validExitNormal )
1015  {
1016  const G4RotationMatrix* rot = motherPhysical->GetRotation();
1017  if ( rot )
1018  {
1019  exitNormal *= rot->inverse();
1020  }
1021  }
1022  }
1023  else
1024  {
1025  validExitNormal = false;
1026  // calculatedExitNormal= false;
1027  }
1028  }
1029 #endif
1030 
1031 
1032  G4bool daughterDeterminedStep=false;
1033  G4ThreeVector daughtNormRepCrd;
1034  // Exit normal of daughter transformed to
1035  // the coordinate system of Replica (i.e. last depth)
1036 
1037  //
1038  // Compute daughter safeties & intersections
1039  //
1040  localNoDaughters = repLogical->GetNoDaughters();
1041  for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1042  {
1043  samplePhysical = repLogical->GetDaughter(sampleNo);
1044  if ( samplePhysical!=blockedExitedVol )
1045  {
1046  G4ThreeVector localExitNorm;
1047  G4ThreeVector normReplicaCoord;
1048 
1049  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
1050  samplePhysical->GetTranslation());
1051  sampleTf.Invert();
1052  const G4ThreeVector samplePoint =
1053  sampleTf.TransformPoint(localPoint);
1054  const G4VSolid* sampleSolid =
1055  samplePhysical->GetLogicalVolume()->GetSolid();
1056  const G4double sampleSafetyDistance =
1057  sampleSolid->DistanceToIn(samplePoint);
1058  if ( sampleSafetyDistance<ourSafety )
1059  {
1060  ourSafety = sampleSafetyDistance;
1061  }
1062  if ( sampleSafetyDistance<=ourStep )
1063  {
1064  sampleDirection = sampleTf.TransformAxis(localDirection);
1065  const G4double sampleStepDistance =
1066  sampleSolid->DistanceToIn(samplePoint,sampleDirection);
1067  if ( sampleStepDistance<=ourStep )
1068  {
1069  daughterDeterminedStep= true;
1070 
1071  ourStep = sampleStepDistance;
1072  entering = true;
1073  exiting = false;
1074  *pBlockedPhysical = samplePhysical;
1075  blockedReplicaNo = sampleNo;
1076 
1077 #ifdef DAUGHTER_NORMAL_ALSO
1078  // This norm can be calculated later, if needed daughter is available
1079  localExitNorm = sampleSolid->SurfaceNormal(samplePoint);
1080  daughtNormRepCrd = sampleTf.Inverse().TransformAxis(localExitNorm);
1081 #endif
1082 
1083 #ifdef G4VERBOSE
1084  // Check to see that the resulting point is indeed in/on volume.
1085  // This check could eventually be made only for successful candidate.
1086 
1087  if ( ( fCheck ) && ( sampleStepDistance < kInfinity ) )
1088  {
1089  G4ThreeVector intersectionPoint;
1090  intersectionPoint= samplePoint
1091  + sampleStepDistance * sampleDirection;
1092  EInside insideIntPt= sampleSolid->Inside(intersectionPoint);
1093  if ( insideIntPt != kSurface )
1094  {
1095  G4int oldcoutPrec = G4cout.precision(16);
1096  std::ostringstream message;
1097  message << "Navigator gets conflicting response from Solid."
1098  << G4endl
1099  << " Inaccurate DistanceToIn for solid "
1100  << sampleSolid->GetName() << G4endl
1101  << " Solid gave DistanceToIn = "
1102  << sampleStepDistance << " yet returns " ;
1103  if ( insideIntPt == kInside )
1104  message << "-kInside-";
1105  else if ( insideIntPt == kOutside )
1106  message << "-kOutside-";
1107  else
1108  message << "-kSurface-";
1109  message << " for this point !" << G4endl
1110  << " Point = " << intersectionPoint << G4endl;
1111  if ( insideIntPt != kInside )
1112  message << " DistanceToIn(p) = "
1113  << sampleSolid->DistanceToIn(intersectionPoint)
1114  << G4endl;
1115  if ( insideIntPt != kOutside )
1116  message << " DistanceToOut(p) = "
1117  << sampleSolid->DistanceToOut(intersectionPoint);
1118  G4Exception("G4ReplicaNavigation::ComputeStep()",
1119  "GeomNav1002", JustWarning, message);
1120  G4cout.precision(oldcoutPrec);
1121  }
1122  }
1123 #endif
1124  }
1125  }
1126  }
1127  }
1128 
1129  calculatedExitNormal &= (!daughterDeterminedStep);
1130 
1131 #ifdef DAUGHTER_NORMAL_ALSO
1132  if( daughterDeterminedStep )
1133  {
1134  // G4ThreeVector daughtNormGlobal =
1135  // GlobalToLastDepth.Inverse().TransformAxis(daughtNormRepCrd);
1136  // ==> Can calculate it, but have no way to transmit it to caller (for now)
1137 
1138  exitNormalVector=globalToLocalTop.Inverse().TransformAxis(daughtNormGlobal);
1139  validExitNormal= false; // Entering daughter - never convex for parent
1140 
1141  calculatedExitNormal= true;
1142  }
1143  // calculatedExitNormal= true; // Force it to true -- dubious
1144 #endif
1145 
1146  newSafety = ourSafety;
1147  return ourStep;
1148 }
1149 
1150 // ********************************************************************
1151 // ComputeSafety
1152 //
1153 // Compute the isotropic distance to current volume's boundaries
1154 // and to daughter volumes.
1155 // ********************************************************************
1156 //
1157 G4double
1159  const G4ThreeVector &localPoint,
1160  G4NavigationHistory &history,
1161  const G4double )
1162 {
1163  G4VPhysicalVolume *repPhysical, *motherPhysical;
1164  G4VPhysicalVolume *samplePhysical, *blockedExitedVol=0;
1165  G4LogicalVolume *repLogical;
1166  G4VSolid *motherSolid;
1167  G4ThreeVector repPoint;
1168  G4double ourSafety=kInfinity;
1169  G4double sampleSafety;
1170  G4int localNoDaughters, sampleNo;
1171  G4int depth;
1172 
1173  repPhysical = history.GetTopVolume();
1174  repLogical = repPhysical->GetLogicalVolume();
1175 
1176  //
1177  // Compute intersection with replica boundaries & replica safety
1178  //
1179 
1180  sampleSafety = DistanceToOut(history.GetTopVolume(),
1181  history.GetTopReplicaNo(),
1182  localPoint);
1183  if ( sampleSafety<ourSafety )
1184  {
1185  ourSafety = sampleSafety;
1186  }
1187 
1188  depth = history.GetDepth()-1;
1189 
1190  // Loop checking, 07.10.2016, J.Apostolakis -- need to add: assert(depth>0)
1191  while ( history.GetVolumeType(depth)==kReplica )
1192  {
1193  repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1194  sampleSafety = DistanceToOut(history.GetVolume(depth),
1195  history.GetReplicaNo(depth),
1196  repPoint);
1197  if ( sampleSafety<ourSafety )
1198  {
1199  ourSafety = sampleSafety;
1200  }
1201  depth--;
1202  }
1203 
1204  // Compute mother safety & intersection
1205  //
1206  repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1207  motherPhysical = history.GetVolume(depth);
1208  motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
1209  sampleSafety = motherSolid->DistanceToOut(repPoint);
1210 
1211  if ( sampleSafety<ourSafety )
1212  {
1213  ourSafety = sampleSafety;
1214  }
1215 
1216  // Compute daughter safeties & intersections
1217  //
1218  localNoDaughters = repLogical->GetNoDaughters();
1219  for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
1220  {
1221  samplePhysical = repLogical->GetDaughter(sampleNo);
1222  if ( samplePhysical!=blockedExitedVol )
1223  {
1224  G4AffineTransform sampleTf(samplePhysical->GetRotation(),
1225  samplePhysical->GetTranslation());
1226  sampleTf.Invert();
1227  const G4ThreeVector samplePoint =
1228  sampleTf.TransformPoint(localPoint);
1229  const G4VSolid *sampleSolid =
1230  samplePhysical->GetLogicalVolume()->GetSolid();
1231  const G4double sampleSafetyDistance =
1232  sampleSolid->DistanceToIn(samplePoint);
1233  if ( sampleSafetyDistance<ourSafety )
1234  {
1235  ourSafety = sampleSafetyDistance;
1236  }
1237  }
1238  }
1239  return ourSafety;
1240 }
1241 
1242 // ********************************************************************
1243 // BackLocate
1244 // ********************************************************************
1245 //
1246 EInside
1248  const G4ThreeVector &globalPoint,
1249  G4ThreeVector &localPoint,
1250  const G4bool &exiting,
1251  G4bool &notKnownInside ) const
1252 {
1253  G4VPhysicalVolume *pNRMother=0;
1254  G4VSolid *motherSolid;
1255  G4ThreeVector repPoint, goodPoint;
1256  G4int mdepth, depth, cdepth;
1257  EInside insideCode;
1258 
1259  cdepth = history.GetDepth();
1260 
1261  // Find non replicated mother
1262  //
1263  for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1264  {
1265  if ( history.GetVolumeType(mdepth)!=kReplica )
1266  {
1267  pNRMother = history.GetVolume(mdepth);
1268  break;
1269  }
1270  }
1271 
1272  if( pNRMother==0 )
1273  {
1274  // All the tree of mother volumes were Replicas.
1275  // This is an error, as the World volume must be a Placement
1276  //
1277  G4Exception("G4ReplicaNavigation::BackLocate()", "GeomNav0002",
1278  FatalException, "The World volume must be a Placement!");
1279  return kInside;
1280  }
1281 
1282  motherSolid = pNRMother->GetLogicalVolume()->GetSolid();
1283  goodPoint = history.GetTransform(mdepth).TransformPoint(globalPoint);
1284  insideCode = motherSolid->Inside(goodPoint);
1285  if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1286  {
1287  // Outside mother -> back up to mother level
1288  // Locate.. in Navigator will back up one more level
1289  // localPoint not required
1290  //
1291  history.BackLevel(cdepth-mdepth);
1292  // localPoint = goodPoint;
1293  }
1294  else
1295  {
1296  notKnownInside = false;
1297 
1298  // Still within replications
1299  // Check down: if on outside stop at this level
1300  //
1301  for ( depth=mdepth+1; depth<cdepth; depth++)
1302  {
1303  repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1304  insideCode = Inside(history.GetVolume(depth),
1305  history.GetReplicaNo(depth),
1306  repPoint);
1307  if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1308  {
1309  localPoint = goodPoint;
1310  history.BackLevel(cdepth-depth);
1311  return insideCode;
1312  }
1313  else
1314  {
1315  goodPoint = repPoint;
1316  }
1317  }
1318  localPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1319  insideCode = Inside(history.GetVolume(depth),
1320  history.GetReplicaNo(depth),
1321  localPoint);
1322  // If outside level, set localPoint = coordinates in reference system
1323  // of *previous* level - location code in navigator will back up one
1324  // level [And also manage blocking]
1325  //
1326  if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1327  {
1328  localPoint = goodPoint;
1329  }
1330  }
1331  return insideCode;
1332 }
Definition: geomdefs.hh:54
G4String GetName() const
G4VPhysicalVolume * GetTopVolume() const
double perp2() const
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
double dot(const Hep3Vector &) const
G4AffineTransform Inverse() const
static const G4double d2
EInside Inside(const G4VPhysicalVolume *pVol, const G4int replicaNo, const G4ThreeVector &localPoint) const
G4int GetDepth() const
G4double GetSurfaceTolerance() const
G4VSolid * GetSolid() const
G4VPhysicalVolume * GetDaughter(const G4int i) const
#define width
const G4RotationMatrix * GetRotation() const
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 setY(double)
double z() const
void setZ(double)
void setX(double)
HepRotation inverse() const
void DumpInfo() const
G4int GetTopReplicaNo() const
G4ThreeVector exitNormal
tuple b
Definition: test.py:12
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
double phi() const
G4int GetNoDaughters() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4ThreeVector & GetTranslation() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4LogicalVolume * GetLogicalVolume() const
const G4AffineTransform & GetTransform(G4int n) const
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
tuple t1
Definition: plottest35.py:33
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
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
tuple c
Definition: test.py:13
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
EInside BackLocate(G4NavigationHistory &history, const G4ThreeVector &globalPoint, G4ThreeVector &localPoint, const G4bool &exiting, G4bool &notKnownInside) const
double perp() const
G4double GetAngularTolerance() const
static G4GeometryTolerance * GetInstance()