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