Geant4  10.00.p02
G4VIntersectionLocator.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 // $Id: G4VIntersectionLocator.cc 66872 2013-01-15 01:25:57Z japost $
27 //
28 // Class G4VIntersectionLocator implementation
29 //
30 // 27.10.08 - John Apostolakis, Tatiana Nikitina.
31 // ---------------------------------------------------------------------------
32 
33 #include <iomanip>
34 #include <sstream>
35 
36 #include "globals.hh"
37 #include "G4ios.hh"
38 #include "G4SystemOfUnits.hh"
40 #include "G4GeometryTolerance.hh"
41 
43 //
44 // Constructor
45 //
47  fUseNormalCorrection(false),
48  fiNavigator( theNavigator ),
49  fiChordFinder( 0 ), // Not set - overridden at each step
50  fiEpsilonStep( -1.0 ), // Out of range - overridden at each step
51  fiDeltaIntersection( -1.0 ), // Out of range - overridden at each step
52  fiUseSafety(false), // Default - overridden at each step
53  fpTouchable(0)
54 {
56  fVerboseLevel = 0;
58 }
59 
61 //
62 // Destructor.
63 //
65 {
66  delete fHelpingNavigator;
67  delete fpTouchable;
68 }
69 
71 //
72 // Dumps status of propagator.
73 //
74 void
76  const G4FieldTrack& CurrentFT,
77  G4double requestStep,
78  G4double safety,
79  G4int stepNo)
80 {
81  const G4int verboseLevel= fVerboseLevel;
82  const G4ThreeVector StartPosition = StartFT.GetPosition();
83  const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
84  const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
85  const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
86 
87  G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
88  G4int oldprc; // cout/cerr precision settings
89 
90  if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
91  {
92  oldprc = G4cout.precision(4);
93  G4cout << std::setw( 6) << " "
94  << std::setw( 25) << " Current Position and Direction" << " "
95  << G4endl;
96  G4cout << std::setw( 5) << "Step#"
97  << std::setw(10) << " s " << " "
98  << std::setw(10) << "X(mm)" << " "
99  << std::setw(10) << "Y(mm)" << " "
100  << std::setw(10) << "Z(mm)" << " "
101  << std::setw( 7) << " N_x " << " "
102  << std::setw( 7) << " N_y " << " "
103  << std::setw( 7) << " N_z " << " " ;
104  G4cout << std::setw( 7) << " Delta|N|" << " "
105  << std::setw( 9) << "StepLen" << " "
106  << std::setw(12) << "StartSafety" << " "
107  << std::setw( 9) << "PhsStep" << " ";
108  G4cout << G4endl;
109  G4cout.precision(oldprc);
110  }
111  if((stepNo == 0) && (verboseLevel <=3))
112  {
113  // Recurse to print the start values
114  //
115  printStatus( StartFT, StartFT, -1.0, safety, -1);
116  }
117  if( verboseLevel <= 3 )
118  {
119  if( stepNo >= 0)
120  {
121  G4cout << std::setw( 4) << stepNo << " ";
122  }
123  else
124  {
125  G4cout << std::setw( 5) << "Start" ;
126  }
127  oldprc = G4cout.precision(8);
128  G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
129  G4cout << std::setw(10) << CurrentPosition.x() << " "
130  << std::setw(10) << CurrentPosition.y() << " "
131  << std::setw(10) << CurrentPosition.z() << " ";
132  G4cout.precision(4);
133  G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
134  << std::setw( 7) << CurrentUnitVelocity.y() << " "
135  << std::setw( 7) << CurrentUnitVelocity.z() << " ";
136  G4cout.precision(3);
137  G4cout << std::setw( 7)
138  << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag()
139  << " ";
140  G4cout << std::setw( 9) << step_len << " ";
141  G4cout << std::setw(12) << safety << " ";
142  if( requestStep != -1.0 )
143  {
144  G4cout << std::setw( 9) << requestStep << " ";
145  }
146  else
147  {
148  G4cout << std::setw( 9) << "Init/NotKnown" << " ";
149  }
150  G4cout << G4endl;
151  G4cout.precision(oldprc);
152  }
153  else // if( verboseLevel > 3 )
154  {
155  // Multi-line output
156 
157  G4cout << "Step taken was " << step_len
158  << " out of PhysicalStep= " << requestStep << G4endl;
159  G4cout << "Final safety is: " << safety << G4endl;
160  G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
161  << G4endl;
162  G4cout << G4endl;
163  }
164 }
165 
167 //
168 // ReEstimateEndPoint.
169 //
171 ReEstimateEndpoint( const G4FieldTrack& CurrentStateA,
172  const G4FieldTrack& EstimatedEndStateB,
173  G4double linearDistSq,
174  G4double curveDist )
175 {
176  G4FieldTrack newEndPoint( CurrentStateA );
178 
179  G4FieldTrack retEndPoint( CurrentStateA );
180  G4bool goodAdvance;
181  G4int itrial=0;
182  const G4int no_trials= 20;
183 
184  G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
185  do
186  {
187  G4double currentCurveLen= newEndPoint.GetCurveLength();
188  G4double advanceLength= endCurveLen - currentCurveLen ;
189  if (std::abs(advanceLength)<kCarTolerance)
190  {
191  goodAdvance=true;
192  }
193  else{
194  goodAdvance=
195  integrDriver->AccurateAdvance(newEndPoint, advanceLength,
197  }
198  }
199  while( !goodAdvance && (++itrial < no_trials) );
200 
201  if( goodAdvance )
202  {
203  retEndPoint= newEndPoint;
204  }
205  else
206  {
207  retEndPoint= EstimatedEndStateB; // Could not improve without major work !!
208  }
209 
210  // All the work is done
211  // below are some diagnostics only -- before the return!
212  //
213  static const G4String MethodName("G4VIntersectionLocator::ReEstimateEndpoint");
214 
215 #ifdef G4VERBOSE
216  G4int latest_good_trials=0;
217  if( itrial > 1)
218  {
219  if( fVerboseLevel > 0 )
220  {
221  G4cout << MethodName << " called - goodAdv= " << goodAdvance
222  << " trials = " << itrial
223  << " previous good= " << latest_good_trials
224  << G4endl;
225  }
226  latest_good_trials=0;
227  }
228  else
229  {
230  latest_good_trials++;
231  }
232 #endif
233 
234 #ifdef G4DEBUG_FIELD
235  G4double lengthDone = newEndPoint.GetCurveLength()
236  - CurrentStateA.GetCurveLength();
237  if( !goodAdvance )
238  {
239  if( fVerboseLevel >= 3 )
240  {
241  G4cout << MethodName << "> AccurateAdvance failed " ;
242  G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
243  G4cout << " It went only " << lengthDone << " instead of " << curveDist
244  << " -- a difference of " << curveDist - lengthDone << G4endl;
245  G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
246  << G4endl;
247  }
248  }
249 
250  static G4int noInaccuracyWarnings = 0;
251  G4int maxNoWarnings = 10;
252  if ( (noInaccuracyWarnings < maxNoWarnings )
253  || (fVerboseLevel > 1) )
254  {
255  G4cerr << "G4PropagatorInField::LocateIntersectionPoint():"
256  << G4endl
257  << " Warning: Integration inaccuracy requires"
258  << " an adjustment in the step's endpoint." << G4endl
259  << " Two mid-points are further apart than their"
260  << " curve length difference" << G4endl
261  << " Dist = " << std::sqrt(linearDistSq)
262  << " curve length = " << curveDist << G4endl;
263  G4cerr << " Correction applied is "
264  << (newEndPoint.GetPosition()-EstimatedEndStateB.GetPosition()).mag()
265  << G4endl;
266  }
267 #else
268  // Statistics on the RMS value of the corrections
269 
270  static G4ThreadLocal G4int noCorrections=0;
271  static G4ThreadLocal G4double sumCorrectionsSq = 0;
272  noCorrections++;
273  if( goodAdvance )
274  {
275  sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
276  newEndPoint.GetPosition()).mag2();
277  }
278  linearDistSq -= curveDist; // To use linearDistSq ... !
279 #endif
280 
281  return retEndPoint;
282 }
283 
285 //
286 // Method for finding SurfaceNormal of Intersecting Solid
287 //
289 GetLocalSurfaceNormal(const G4ThreeVector& CurrentE_Point, G4bool& validNormal)
290 {
291  G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0.0));
292  G4VPhysicalVolume* located;
293 
294  validNormal = false;
295  fHelpingNavigator->SetWorldVolume(GetNavigatorFor()->GetWorldVolume());
296  located = fHelpingNavigator->LocateGlobalPointAndSetup( CurrentE_Point );
297 
298  delete fpTouchable;
300 
301  // To check if we can use GetGlobalExitNormal()
302  //
303  G4ThreeVector localPosition = fpTouchable->GetHistory()
304  ->GetTopTransform().TransformPoint(CurrentE_Point);
305 
306  // Issue: in the case of coincident surfaces, this version does not recognise
307  // which side you are located onto (can return vector with wrong sign.)
308  // TO-DO: use direction (of chord) to identify volume we will be "entering"
309 
310  if( located != 0)
311  {
312  G4LogicalVolume* pLogical= located->GetLogicalVolume();
313  G4VSolid* pSolid;
314 
315  if( (pLogical != 0) && ( (pSolid=pLogical->GetSolid()) !=0 ) )
316  {
317  // G4bool goodPoint, nearbyPoint;
318  // G4int numGoodPoints, numNearbyPoints; // --> use for stats
319  if ( ( pSolid->Inside(localPosition)==kSurface )
320  || ( pSolid->DistanceToOut(localPosition) < 1000.0 * kCarTolerance )
321  )
322  {
323  Normal = pSolid->SurfaceNormal(localPosition);
324  validNormal = true;
325 
326 #ifdef G4DEBUG_FIELD
327  if( std::fabs(Normal.mag2() - 1.0 ) > perMille)
328  {
329  G4cerr << "PROBLEM in G4VIntersectionLocator::GetLocalSurfaceNormal."
330  << G4endl;
331  G4cerr << " Normal is not unit - mag=" << Normal.mag() << G4endl;
332  G4cerr << " at trial local point " << CurrentE_Point << G4endl;
333  G4cerr << " Solid is " << *pSolid << G4endl;
334  }
335 #endif
336  }
337  }
338  }
339 
340  return Normal;
341 }
342 
344 //
345 // Adjustment of Found Intersection
346 //
349  const G4ThreeVector& CurrentE_Point,
350  const G4ThreeVector& CurrentF_Point,
351  const G4ThreeVector& MomentumDir,
352  const G4bool IntersectAF,
353  G4ThreeVector& IntersectionPoint, // I/O
354  G4double& NewSafety, // I/O
355  G4double& fPreviousSafety, // I/O
356  G4ThreeVector& fPreviousSftOrigin )// I/O
357 {
358  G4double dist,lambda;
359  G4ThreeVector Normal, NewPoint, Point_G;
360  G4bool goodAdjust=false, Intersects_FP=false, validNormal=false;
361 
362  // Get SurfaceNormal of Intersecting Solid
363  //
364  Normal = GetGlobalSurfaceNormal(CurrentE_Point,validNormal);
365  if(!validNormal) { return false; }
366 
367  // Intersection between Line and Plane
368  //
369  G4double n_d_m = Normal.dot(MomentumDir);
370  if ( std::abs(n_d_m)>kCarTolerance )
371  {
372 #ifdef G4VERBOSE
373  if ( fVerboseLevel>1 )
374  {
375  G4cerr << "WARNING - "
376  << "G4VIntersectionLocator::AdjustementOfFoundIntersection()"
377  << G4endl
378  << " No intersection. Parallels lines!" << G4endl;
379  }
380 #endif
381  lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
382 
383  // New candidate for Intersection
384  //
385  NewPoint = CurrentF_Point+lambda*MomentumDir;
386 
387  // Distance from CurrentF to Calculated Intersection
388  //
389  dist = std::abs(lambda);
390 
391  if ( dist<kCarTolerance*0.001 ) { return false; }
392 
393  // Calculation of new intersection point on the path.
394  //
395  if ( IntersectAF ) // First part intersects
396  {
397  G4double stepLengthFP;
398  G4ThreeVector Point_P = CurrentA_Point;
400  Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
401  fPreviousSafety, fPreviousSftOrigin,
402  stepLengthFP, Point_G );
403 
404  }
405  else // Second part intersects
406  {
407  G4double stepLengthFP;
409  Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
410  fPreviousSafety, fPreviousSftOrigin,
411  stepLengthFP, Point_G );
412  }
413  if ( Intersects_FP )
414  {
415  goodAdjust = true;
416  IntersectionPoint = Point_G;
417  }
418  }
419 
420  return goodAdjust;
421 }
422 
425  G4bool& validNormal) // const
426 {
427  G4ThreeVector NormalAtEntry; // ( -10. , -10., -10. );
428 
429  G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
430  G4bool validNormalLast;
431 
432  // Relies on a call to Navigator::ComputeStep in IntersectChord before
433  // this call
434  //
435  NormalAtEntryLast = GetLastSurfaceNormal( CurrentInt_Point, validNormalLast );
436  // May return valid=false in cases, including
437  // - if the candidate volume was not found (eg exiting world), or
438  // - a replica was involved -- determined the step size.
439  // (This list is not complete.)
440 
441 #ifdef G4DEBUG_FIELD
442  if ( validNormalLast
443  && ( std::fabs(NormalAtEntryLast.mag2() - 1.0) > perThousand ) )
444  {
445  std::ostringstream message;
446  message << "G4VIntersectionLocator::GetSurfaceNormal -- identified problem."
447  << G4endl;
448  message << "PROBLEM: Normal is not unit - magnitude = "
449  << NormalAtEntryLast.mag() << G4endl;
450  message << " at trial intersection point " << CurrentInt_Point << G4endl;
451  message << " Obtained from Get *Last* Surface Normal." << G4endl;
452  G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
453  "GeomNav1002", JustWarning, message);
454  }
455 #endif
456 
457  if( validNormalLast )
458  {
459  NormalAtEntry=NormalAtEntryLast;
460  }
461  validNormal = validNormalLast;
462 
463  return NormalAtEntry;
464 }
465 
467 GetGlobalSurfaceNormal(const G4ThreeVector& CurrentE_Point,
468  G4bool& validNormal)
469 {
470  G4ThreeVector localNormal=
471  GetLocalSurfaceNormal( CurrentE_Point, validNormal );
472  G4AffineTransform localToGlobal= // Must use the same Navigator !!
474  G4ThreeVector globalNormal =
475  localToGlobal.TransformAxis( localNormal );
476 
477 #ifdef G4DEBUG_FIELD
478  if( validNormal && ( std::fabs(globalNormal.mag2() - 1.0) > perThousand ) )
479  {
480  std::ostringstream message;
481  message << "**************************************************************"
482  << G4endl;
483  message << " Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
484  << G4endl;
485  message << " * Constituents: " << G4endl;
486  message << " Local Normal= " << localNormal << G4endl;
487  message << " Transform: " << G4endl
488  << " Net Translation= " << localToGlobal.NetTranslation()
489  << G4endl
490  << " Net Rotation = " << localToGlobal.NetRotation()
491  << G4endl;
492  message << " * Result: " << G4endl;
493  message << " Global Normal= " << localNormal << G4endl;
494  message << "**************************************************************"
495  << G4endl;
496  G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
497  "GeomNav1002", JustWarning, message);
498  }
499 #endif
500 
501  return globalNormal;
502 }
503 
506  G4bool& normalIsValid) const
507 {
508  G4ThreeVector normalVec;
509  G4bool validNorm;
510  normalVec = fiNavigator->GetGlobalExitNormal( intersectPoint, &validNorm );
511  normalIsValid= validNorm;
512 
513  return normalVec;
514 }
515 
517  const G4ThreeVector& ChordAB_v,
518  const G4ThreeVector& ChordEF_v,
519  const G4ThreeVector& NewMomentumDir,
520  const G4ThreeVector& NormalAtEntry,
521  G4bool validNormal )
522 {
523  G4double ABchord_length = ChordAB_v.mag();
524  G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ) ;
525  G4double MomDir_dot_ABchord;
526  MomDir_dot_ABchord= (1.0 / ABchord_length) * NewMomentumDir.dot( ChordAB_v );
527 
528  std::ostringstream outStream;
529  outStream // G4cout
530  << std::setw(6) << " Step# "
531  << std::setw(17) << " |ChordEF|(mag)" << " "
532  << std::setw(18) << " uMomentum.Normal" << " "
533  << std::setw(18) << " uMomentum.ABdir " << " "
534  << std::setw(16) << " AB-dist " << " "
535  << " Chord Vector (EF) "
536  << G4endl;
537  outStream.precision(7);
538  outStream // G4cout
539  << " " << std::setw(5) << step_no
540  << " " << std::setw(18) << ChordEF_v.mag()
541  << " " << std::setw(18) << MomDir_dot_Norm
542  << " " << std::setw(18) << MomDir_dot_ABchord
543  << " " << std::setw(12) << ABchord_length
544  << " " << ChordEF_v
545  << G4endl;
546  outStream // G4cout
547  << " MomentumDir= " << " " << NewMomentumDir
548  << " Normal at Entry E= " << NormalAtEntry
549  << " AB chord = " << ChordAB_v
550  << G4endl;
551  G4cout << outStream.str(); // ostr_verbose;
552 
553  if( ( std::fabs(NormalAtEntry.mag2() - 1.0) > perThousand ) )
554  {
555  G4cerr << " PROBLEM in G4VIntersectionLocator::ReportTrialStep " << G4endl
556  << " Normal is not unit - mag=" << NormalAtEntry.mag()
557  << " ValidNormalAtE = " << validNormal
558  << G4endl;
559  }
560  return;
561 }
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetGlobalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4double GetSurfaceTolerance() const
const G4ThreeVector & GetMomentumDir() const
G4ThreeVector GetLastSurfaceNormal(const G4ThreeVector &intersectPoint, G4bool &validNormal) const
G4ThreeVector NetTranslation() const
static const double perThousand
Definition: G4SIunits.hh:297
const G4AffineTransform GetLocalToGlobalTransform() const
#define G4ThreadLocal
Definition: tls.hh:52
int G4int
Definition: G4Types.hh:78
G4double GetEpsilonStepFor()
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
virtual EInside Inside(const G4ThreeVector &p) const =0
bool G4bool
Definition: G4Types.hh:79
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4Navigator * GetNavigatorFor()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4RotationMatrix NetRotation() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4TouchableHistory * CreateTouchableHistory() const
G4LogicalVolume * GetLogicalVolume() const
virtual G4ThreeVector GetGlobalExitNormal(const G4ThreeVector &point, G4bool *valid)
G4VIntersectionLocator(G4Navigator *theNavigator)
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step)
void SetWorldVolume(G4VPhysicalVolume *pWorld)
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:118
G4ThreeVector GetLocalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
#define G4endl
Definition: G4ios.hh:61
const G4AffineTransform & GetTopTransform() const
const G4NavigationHistory * GetHistory() const
G4ThreeVector GetMomentum() const
double G4double
Definition: G4Types.hh:76
G4TouchableHistory * fpTouchable
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4MagInt_Driver * GetIntegrationDriver()
G4ChordFinder * GetChordFinderFor()
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
static G4GeometryTolerance * GetInstance()
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:550
G4VSolid * GetSolid() const
G4bool AdjustmentOfFoundIntersection(const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
G4GLOB_DLL std::ostream G4cerr