ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VIntersectionLocator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VIntersectionLocator.cc
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 // Class G4VIntersectionLocator implementation
27 //
28 // 27.10.08 - John Apostolakis, Tatiana Nikitina.
29 // ---------------------------------------------------------------------------
30 
31 #include <iomanip>
32 #include <sstream>
33 
34 #include "globals.hh"
35 #include "G4ios.hh"
36 #include "G4AutoDelete.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4VIntersectionLocator.hh"
39 #include "G4GeometryTolerance.hh"
40 
42 //
43 // Constructor
44 //
45 G4VIntersectionLocator::G4VIntersectionLocator(G4Navigator* theNavigator)
46  : fiNavigator(theNavigator)
47 {
48  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
49 
50  if( fiNavigator->GetExternalNavigation() == nullptr )
51  {
52  fHelpingNavigator = new G4Navigator();
53  }
54  else // Must clone the navigator, together with External Navigation
55  {
56  fHelpingNavigator = fiNavigator->Clone();
57  }
58 }
59 
61 //
62 // Destructor.
63 //
64 G4VIntersectionLocator::~G4VIntersectionLocator()
65 {
66  delete fHelpingNavigator;
67  delete fpTouchable;
68 }
69 
71 //
72 // Dump status of propagator to cout (old method)
73 //
74 void
75 G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
76  const G4FieldTrack& CurrentFT,
77  G4double requestStep,
78  G4double safety,
79  G4int stepNo)
80 {
81  std::ostringstream os;
82  printStatus( StartFT,CurrentFT,requestStep,safety,stepNo,os,fVerboseLevel);
83  G4cout << os.str();
84 }
85 
87 //
88 // Dumps status of propagator.
89 //
90 void
91 G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
92  const G4FieldTrack& CurrentFT,
93  G4double requestStep,
94  G4double safety,
95  G4int stepNo,
96  std::ostream& os,
97  G4int verboseLevel)
98 {
99  // const G4int verboseLevel= fVerboseLevel;
100  const G4ThreeVector StartPosition = StartFT.GetPosition();
101  const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
102  const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
103  const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
104 
105  G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
106  G4int oldprc; // cout/cerr precision settings
107 
108  if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
109  {
110  oldprc = os.precision(4);
111  os << std::setw( 6) << " "
112  << std::setw( 25) << " Current Position and Direction" << " "
113  << G4endl;
114  os << std::setw( 5) << "Step#"
115  << std::setw(10) << " s " << " "
116  << std::setw(10) << "X(mm)" << " "
117  << std::setw(10) << "Y(mm)" << " "
118  << std::setw(10) << "Z(mm)" << " "
119  << std::setw( 7) << " N_x " << " "
120  << std::setw( 7) << " N_y " << " "
121  << std::setw( 7) << " N_z " << " " ;
122  os << std::setw( 7) << " Delta|N|" << " "
123  << std::setw( 9) << "StepLen" << " "
124  << std::setw(12) << "StartSafety" << " "
125  << std::setw( 9) << "PhsStep" << " ";
126  os << G4endl;
127  os.precision(oldprc);
128  }
129  if((stepNo == 0) && (verboseLevel <=3))
130  {
131  // Recurse to print the start values
132  //
133  printStatus( StartFT, StartFT, -1.0, safety, -1, os, verboseLevel);
134  }
135  if( verboseLevel <= 3 )
136  {
137  if( stepNo >= 0)
138  {
139  os << std::setw( 4) << stepNo << " ";
140  }
141  else
142  {
143  os << std::setw( 5) << "Start" ;
144  }
145  oldprc = os.precision(8);
146  os << std::setw(10) << CurrentFT.GetCurveLength() << " ";
147  os << std::setw(10) << CurrentPosition.x() << " "
148  << std::setw(10) << CurrentPosition.y() << " "
149  << std::setw(10) << CurrentPosition.z() << " ";
150  os.precision(4);
151  os << std::setw( 7) << CurrentUnitVelocity.x() << " "
152  << std::setw( 7) << CurrentUnitVelocity.y() << " "
153  << std::setw( 7) << CurrentUnitVelocity.z() << " ";
154  os.precision(3);
155  os << std::setw( 7)
156  << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag()
157  << " ";
158  os << std::setw( 9) << step_len << " ";
159  os << std::setw(12) << safety << " ";
160  if( requestStep != -1.0 )
161  {
162  os << std::setw( 9) << requestStep << " ";
163  }
164  else
165  {
166  os << std::setw( 9) << "Init/NotKnown" << " ";
167  }
168  os << G4endl;
169  os.precision(oldprc);
170  }
171  else // if( verboseLevel > 3 )
172  {
173  // Multi-line output
174 
175  os << "Step taken was " << step_len
176  << " out of PhysicalStep= " << requestStep << G4endl;
177  os << "Final safety is: " << safety << G4endl;
178  os << "Chord length = " << (CurrentPosition-StartPosition).mag()
179  << G4endl;
180  os << G4endl;
181  }
182 }
183 
185 //
186 // ReEstimateEndPoint.
187 //
188 G4FieldTrack G4VIntersectionLocator::
189 ReEstimateEndpoint( const G4FieldTrack& CurrentStateA,
190  const G4FieldTrack& EstimatedEndStateB,
191  G4double , // linearDistSq, // NOT used
192  G4double
193 #ifdef G4DEBUG_FIELD
194  curveDist
195 #endif
196  )
197 {
198  G4FieldTrack newEndPoint( CurrentStateA );
199  auto integrDriver = GetChordFinderFor()->GetIntegrationDriver();
200 
201  G4FieldTrack retEndPoint( CurrentStateA );
202  G4bool goodAdvance;
203  G4int itrial = 0;
204  const G4int no_trials = 20;
205 
206 
207  G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
208 
209  do // Loop checking, 07.10.2016, JA
210  {
211  G4double currentCurveLen = newEndPoint.GetCurveLength();
212  G4double advanceLength = endCurveLen - currentCurveLen ;
213  if (std::abs(advanceLength)<kCarTolerance)
214  {
215  goodAdvance=true;
216  }
217  else
218  {
219  goodAdvance = integrDriver->AccurateAdvance(newEndPoint, advanceLength,
220  GetEpsilonStepFor());
221  }
222  }
223  while( !goodAdvance && (++itrial < no_trials) );
224 
225  if( goodAdvance )
226  {
227  retEndPoint = newEndPoint;
228  }
229  else
230  {
231  retEndPoint = EstimatedEndStateB; // Could not improve without major work !!
232  }
233 
234  // All the work is done
235  // below are some diagnostics only -- before the return!
236  //
237  const G4String MethodName("G4VIntersectionLocator::ReEstimateEndpoint()");
238 
239 #ifdef G4VERBOSE
240  G4int latest_good_trials = 0;
241  if( itrial > 1)
242  {
243  if( fVerboseLevel > 0 )
244  {
245  G4cout << MethodName << " called - goodAdv= " << goodAdvance
246  << " trials = " << itrial
247  << " previous good= " << latest_good_trials
248  << G4endl;
249  }
250  latest_good_trials = 0;
251  }
252  else
253  {
254  ++latest_good_trials;
255  }
256 #endif
257 
258 #ifdef G4DEBUG_FIELD
259  G4double lengthDone = newEndPoint.GetCurveLength()
260  - CurrentStateA.GetCurveLength();
261  if( !goodAdvance )
262  {
263  if( fVerboseLevel >= 3 )
264  {
265  G4cout << MethodName << "> AccurateAdvance failed " ;
266  G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
267  G4cout << " It went only " << lengthDone << " instead of " << curveDist
268  << " -- a difference of " << curveDist - lengthDone << G4endl;
269  G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
270  << G4endl;
271  }
272  }
273  G4double linearDist = ( EstimatedEndStateB.GetPosition()
274  - CurrentStateA.GetPosition() ).mag();
275  static G4int noInaccuracyWarnings = 0;
276  G4int maxNoWarnings = 10;
277  if ( (noInaccuracyWarnings < maxNoWarnings )
278  || (fVerboseLevel > 1) )
279  {
280  G4ThreeVector move = newEndPoint.GetPosition()
281  - EstimatedEndStateB.GetPosition();
282  std::ostringstream message;
283  message.precision(12);
284  message << " Integration inaccuracy requires"
285  << " an adjustment in the step's endpoint." << G4endl
286  << " Two mid-points are further apart than their"
287  << " curve length difference" << G4endl
288  << " Dist = " << linearDist
289  << " curve length = " << curveDist << G4endl;
290  message << " Correction applied is " << move.mag() << G4endl
291  << " Old Estimated B position= "
292  << EstimatedEndStateB.GetPosition() << G4endl
293  << " Recalculated Position= "
294  << newEndPoint.GetPosition() << G4endl
295  << " Change ( new - old ) = " << move;
296  G4Exception("G4VIntersectionLocator::ReEstimateEndpoint()",
297  "GeomNav1002", JustWarning, message);
298  }
299 #else
300  // Statistics on the RMS value of the corrections
301 
302  static G4ThreadLocal G4int noCorrections = 0;
303  static G4ThreadLocal G4double sumCorrectionsSq = 0;
304  ++noCorrections;
305  if( goodAdvance )
306  {
307  sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
308  newEndPoint.GetPosition()).mag2();
309  }
310 #endif
311 
312  return retEndPoint;
313 }
314 
315 
317 //
318 // ReEstimateEndPoint.
319 //
320 // The following values are returned in curveError
321 // 0: Normal - no problem
322 // 1: Unexpected co-incidence - milder mixup
323 // 2: Real mixup - EndB is NOT past StartA
324 // ( ie. StartA's curve-lengh is bigger than EndB's)
325 
326 
327 G4bool G4VIntersectionLocator::
328 CheckAndReEstimateEndpoint( const G4FieldTrack& CurrentStartA,
329  const G4FieldTrack& EstimatedEndB,
330  G4FieldTrack& RevisedEndPoint,
331  G4int& curveError)
332 {
333  G4double linDistSq, curveDist;
334 
335  G4bool recalculated = false;
336  curveError= 0;
337 
338  linDistSq = ( EstimatedEndB.GetPosition()
339  - CurrentStartA.GetPosition() ).mag2();
340  curveDist = EstimatedEndB.GetCurveLength()
341  - CurrentStartA.GetCurveLength();
342  if( (curveDist>=0.0)
343  && (curveDist*curveDist *(1.0+2.0*fiEpsilonStep ) < linDistSq ) )
344  {
345  G4FieldTrack newEndPointFT = EstimatedEndB; // Unused
346 
347  if (curveDist>0.0)
348  {
349  // Re-integrate to obtain a new B
350  RevisedEndPoint = ReEstimateEndpoint( CurrentStartA,
351  EstimatedEndB,
352  linDistSq,
353  curveDist );
354  recalculated = true;
355  }
356  else
357  {
358  // Zero length -> no advance!
359  newEndPointFT = CurrentStartA;
360  recalculated = true;
361  curveError = 1; // Unexpected co-incidence - milder mixup
362 
363  G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
364  "GeomNav1002", JustWarning,
365  "A & B are at equal distance in 2nd half. A & B will coincide." );
366  }
367  }
368 
369  // Sanity check
370  //
371  if( curveDist < 0.0 )
372  {
373  curveError = 2; // Real mixup
374  }
375  return recalculated;
376 }
377 
379 //
380 // Method for finding SurfaceNormal of Intersecting Solid
381 //
382 G4ThreeVector G4VIntersectionLocator::
383 GetLocalSurfaceNormal(const G4ThreeVector& CurrentE_Point, G4bool& validNormal)
384 {
385  G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0.0));
386  G4VPhysicalVolume* located;
387 
388  validNormal = false;
389  fHelpingNavigator->SetWorldVolume(GetNavigatorFor()->GetWorldVolume());
390  located = fHelpingNavigator->LocateGlobalPointAndSetup( CurrentE_Point );
391 
392  delete fpTouchable;
393  fpTouchable = fHelpingNavigator->CreateTouchableHistory();
394 
395  // To check if we can use GetGlobalExitNormal()
396  //
397  G4ThreeVector localPosition = fpTouchable->GetHistory()
398  ->GetTopTransform().TransformPoint(CurrentE_Point);
399 
400  // Issue: in the case of coincident surfaces, this version does not recognise
401  // which side you are located onto (can return vector with wrong sign.)
402  // TO-DO: use direction (of chord) to identify volume we will be "entering"
403 
404  if( located != 0)
405  {
406  G4LogicalVolume* pLogical= located->GetLogicalVolume();
407  G4VSolid* pSolid;
408 
409  if( (pLogical != nullptr) && ( (pSolid=pLogical->GetSolid()) != nullptr ) )
410  {
411  if ( ( pSolid->Inside(localPosition)==kSurface )
412  || ( pSolid->DistanceToOut(localPosition) < 1000.0 * kCarTolerance ) )
413  {
414  Normal = pSolid->SurfaceNormal(localPosition);
415  validNormal = true;
416 
417 #ifdef G4DEBUG_FIELD
418  if( std::fabs(Normal.mag2() - 1.0 ) > CLHEP::perThousand)
419  {
420  G4cerr << "PROBLEM in G4VIntersectionLocator::GetLocalSurfaceNormal."
421  << G4endl;
422  G4cerr << " Normal is not unit - mag=" << Normal.mag() << G4endl;
423  G4cerr << " at trial local point " << CurrentE_Point << G4endl;
424  G4cerr << " Solid is " << *pSolid << G4endl;
425  }
426 #endif
427  }
428  }
429  }
430  return Normal;
431 }
432 
434 //
435 // Adjustment of Found Intersection
436 //
437 G4bool G4VIntersectionLocator::
438 AdjustmentOfFoundIntersection( const G4ThreeVector& CurrentA_Point,
439  const G4ThreeVector& CurrentE_Point,
440  const G4ThreeVector& CurrentF_Point,
441  const G4ThreeVector& MomentumDir,
442  const G4bool IntersectAF,
443  G4ThreeVector& IntersectionPoint, // I/O
444  G4double& NewSafety, // I/O
445  G4double& fPreviousSafety, // I/O
446  G4ThreeVector& fPreviousSftOrigin )// I/O
447 {
448  G4double dist,lambda;
449  G4ThreeVector Normal, NewPoint, Point_G;
450  G4bool goodAdjust = false, Intersects_FP = false, validNormal = false;
451 
452  // Get SurfaceNormal of Intersecting Solid
453  //
454  Normal = GetGlobalSurfaceNormal(CurrentE_Point,validNormal);
455  if(!validNormal) { return false; }
456 
457  // Intersection between Line and Plane
458  //
459  G4double n_d_m = Normal.dot(MomentumDir);
460  if ( std::abs(n_d_m)>kCarTolerance )
461  {
462 #ifdef G4VERBOSE
463  if ( fVerboseLevel>1 )
464  {
465  G4Exception("G4VIntersectionLocator::AdjustmentOfFoundIntersection()",
466  "GeomNav0003", JustWarning,
467  "No intersection. Parallels lines!");
468  }
469 #endif
470  lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
471 
472  // New candidate for Intersection
473  //
474  NewPoint = CurrentF_Point+lambda*MomentumDir;
475 
476  // Distance from CurrentF to Calculated Intersection
477  //
478  dist = std::abs(lambda);
479 
480  if ( dist<kCarTolerance*0.001 ) { return false; }
481 
482  // Calculation of new intersection point on the path.
483  //
484  if ( IntersectAF ) // First part intersects
485  {
486  G4double stepLengthFP;
487  G4ThreeVector Point_P = CurrentA_Point;
488  GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_P);
489  Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
490  fPreviousSafety, fPreviousSftOrigin,
491  stepLengthFP, Point_G );
492 
493  }
494  else // Second part intersects
495  {
496  G4double stepLengthFP;
497  GetNavigatorFor()->LocateGlobalPointWithinVolume(CurrentF_Point );
498  Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
499  fPreviousSafety, fPreviousSftOrigin,
500  stepLengthFP, Point_G );
501  }
502  if ( Intersects_FP )
503  {
504  goodAdjust = true;
505  IntersectionPoint = Point_G;
506  }
507  }
508 
509  return goodAdjust;
510 }
511 
513 //
514 // GetSurfaceNormal.
515 //
516 G4ThreeVector G4VIntersectionLocator::
517 GetSurfaceNormal(const G4ThreeVector& CurrentInt_Point,
518  G4bool& validNormal)
519 {
520  G4ThreeVector NormalAtEntry; // ( -10. , -10., -10. );
521 
522  G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
523  G4bool validNormalLast;
524 
525  // Relies on a call to Navigator::ComputeStep in IntersectChord before
526  // this call
527  //
528  NormalAtEntryLast = GetLastSurfaceNormal( CurrentInt_Point, validNormalLast );
529  // May return valid=false in cases, including
530  // - if the candidate volume was not found (eg exiting world), or
531  // - a replica was involved -- determined the step size.
532  // (This list is not complete.)
533 
534 #ifdef G4DEBUG_FIELD
535  if ( validNormalLast
536  && ( std::fabs(NormalAtEntryLast.mag2() - 1.0) > perThousand ) )
537  {
538  std::ostringstream message;
539  message << "PROBLEM: Normal is not unit - magnitude = "
540  << NormalAtEntryLast.mag() << G4endl;
541  message << " at trial intersection point " << CurrentInt_Point << G4endl;
542  message << " Obtained from Get *Last* Surface Normal.";
543  G4Exception("G4VIntersectionLocator::GetSurfaceNormal()",
544  "GeomNav1002", JustWarning, message);
545  }
546 #endif
547 
548  if( validNormalLast )
549  {
550  NormalAtEntry = NormalAtEntryLast;
551  }
552  validNormal = validNormalLast;
553 
554  return NormalAtEntry;
555 }
556 
558 //
559 // GetGlobalSurfaceNormal.
560 //
561 G4ThreeVector G4VIntersectionLocator::
562 GetGlobalSurfaceNormal(const G4ThreeVector& CurrentE_Point,
563  G4bool& validNormal)
564 {
565  G4ThreeVector localNormal = GetLocalSurfaceNormal(CurrentE_Point,validNormal);
566  G4AffineTransform localToGlobal = // Must use the same Navigator !!
567  fHelpingNavigator->GetLocalToGlobalTransform();
568  G4ThreeVector globalNormal = localToGlobal.TransformAxis( localNormal );
569 
570 #ifdef G4DEBUG_FIELD
571  if( validNormal && ( std::fabs(globalNormal.mag2() - 1.0) > perThousand ) )
572  {
573  std::ostringstream message;
574  message << "**************************************************************"
575  << G4endl;
576  message << " Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
577  << G4endl;
578  message << " * Constituents: " << G4endl;
579  message << " Local Normal= " << localNormal << G4endl;
580  message << " Transform: " << G4endl
581  << " Net Translation= " << localToGlobal.NetTranslation()
582  << G4endl
583  << " Net Rotation = " << localToGlobal.NetRotation()
584  << G4endl;
585  message << " * Result: " << G4endl;
586  message << " Global Normal= " << localNormal << G4endl;
587  message << "**************************************************************";
588  G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
589  "GeomNav1002", JustWarning, message);
590  }
591 #endif
592 
593  return globalNormal;
594 }
595 
597 //
598 // GetLastSurfaceNormal.
599 //
600 G4ThreeVector G4VIntersectionLocator::
601 GetLastSurfaceNormal( const G4ThreeVector& intersectPoint,
602  G4bool& normalIsValid) const
603 {
604  G4ThreeVector normalVec;
605  G4bool validNorm;
606  normalVec = fiNavigator->GetGlobalExitNormal( intersectPoint, &validNorm );
607  normalIsValid = validNorm;
608 
609  return normalVec;
610 }
611 
613 //
614 // ReportTrialStep.
615 //
616 void G4VIntersectionLocator::ReportTrialStep( G4int step_no,
617  const G4ThreeVector& ChordAB_v,
618  const G4ThreeVector& ChordEF_v,
619  const G4ThreeVector& NewMomentumDir,
620  const G4ThreeVector& NormalAtEntry,
621  G4bool validNormal )
622 {
623  G4double ABchord_length = ChordAB_v.mag();
624  G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry );
625  G4double MomDir_dot_ABchord;
626  MomDir_dot_ABchord = (1.0 / ABchord_length) * NewMomentumDir.dot( ChordAB_v );
627 
628  std::ostringstream outStream;
629  outStream << std::setw(6) << " Step# "
630  << std::setw(17) << " |ChordEF|(mag)" << " "
631  << std::setw(18) << " uMomentum.Normal" << " "
632  << std::setw(18) << " uMomentum.ABdir " << " "
633  << std::setw(16) << " AB-dist " << " "
634  << " Chord Vector (EF) "
635  << G4endl;
636  outStream.precision(7);
637  outStream << " " << std::setw(5) << step_no
638  << " " << std::setw(18) << ChordEF_v.mag()
639  << " " << std::setw(18) << MomDir_dot_Norm
640  << " " << std::setw(18) << MomDir_dot_ABchord
641  << " " << std::setw(12) << ABchord_length
642  << " " << ChordEF_v
643  << G4endl;
644  outStream << " MomentumDir= " << " " << NewMomentumDir
645  << " Normal at Entry E= " << NormalAtEntry
646  << " AB chord = " << ChordAB_v
647  << G4endl;
648  G4cout << outStream.str();
649 
650  if( ( std::fabs(NormalAtEntry.mag2() - 1.0) > perThousand ) )
651  {
652  std::ostringstream message;
653  message << "Normal is not unit - mag= " << NormalAtEntry.mag() << G4endl
654  << " ValidNormalAtE = " << validNormal;
655  G4Exception("G4VIntersectionLocator::ReportTrialStep()",
656  "GeomNav1002", JustWarning, message);
657  }
658  return;
659 }
660 
662 //
663 // LocateGlobalPointWithinVolumeAndCheck.
664 //
665 // Locate point using navigator: updates state of Navigator
666 // By default, it assumes that the point is inside the current volume,
667 // and returns true.
668 // In check mode, checks whether the point is *inside* the volume.
669 // If it is inside, it returns true
670 // If not, issues a warning and returns false.
671 //
672 G4bool G4VIntersectionLocator::
673 LocateGlobalPointWithinVolumeAndCheck( const G4ThreeVector& position )
674 {
675  G4bool good = true;
676  G4Navigator* nav = GetNavigatorFor();
677  const G4String
678  MethodName("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()");
679 
680  if( fCheckMode )
681  {
682  G4bool navCheck= nav->IsCheckModeActive(); // Recover original value
683  nav->CheckMode(true);
684 
685  // Identify the current volume
686 
687  G4TouchableHistoryHandle startTH= nav->CreateTouchableHistoryHandle();
688  G4VPhysicalVolume* motherPhys = startTH->GetVolume();
689  G4VSolid* motherSolid = startTH->GetSolid();
690  G4AffineTransform transform = nav->GetGlobalToLocalTransform();
691  G4int motherCopyNo = motherPhys->GetCopyNo();
692 
693  // Let's check that the point is inside the current solid
694  G4ThreeVector localPosition = transform.TransformPoint(position);
695  EInside inMother = motherSolid->Inside( localPosition );
696  if( inMother != kInside )
697  {
698  std::ostringstream message;
699  message << "Position located "
700  << ( inMother == kSurface ? " on Surface " : " outside " )
701  << "expected volume" << G4endl
702  << " Safety (from Outside) = "
703  << motherSolid->DistanceToIn(localPosition);
704  G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
705  "GeomNav1002", JustWarning, message);
706  }
707 
708  // 1. Simple next step - quick relocation and check result.
709  // nav->LocateGlobalPointWithinVolume( position );
710 
711  // 2. Full relocation - to cross-check answer !
712  G4VPhysicalVolume* nextPhysical= nav->LocateGlobalPointAndSetup(position);
713  if( (nextPhysical != motherPhys)
714  || (nextPhysical->GetCopyNo() != motherCopyNo )
715  )
716  {
717  G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
718  "GeomNav1002", JustWarning,
719  "Position located outside expected volume.");
720  }
721  nav->CheckMode(navCheck); // Recover original value
722  }
723  else
724  {
725  nav->LocateGlobalPointWithinVolume( position );
726  }
727  return good;
728 }
729 
731 //
732 // LocateGlobalPointWithinVolumeCheckAndReport.
733 //
734 void G4VIntersectionLocator::
735 LocateGlobalPointWithinVolumeCheckAndReport( const G4ThreeVector& position,
736  const G4String& CodeLocationInfo,
737  G4int /* CheckMode */)
738 {
739  // Save value of Check mode first
740  G4bool oldCheck = GetCheckMode();
741 
742  G4bool ok = LocateGlobalPointWithinVolumeAndCheck( position );
743  if( !ok )
744  {
745  std::ostringstream message;
746  message << "Failed point location." << G4endl
747  << " Code Location info: " << CodeLocationInfo;
748  G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport()",
749  "GeomNav1002", JustWarning, message);
750  }
751 
752  SetCheckMode( oldCheck );
753 }
754 
756 //
757 // ReportReversedPoints.
758 //
759 void G4VIntersectionLocator::
760 ReportReversedPoints( std::ostringstream& msg,
761  const G4FieldTrack& StartPointVel,
762  const G4FieldTrack& EndPointVel,
763  G4double NewSafety, G4double epsStep,
764  const G4FieldTrack& A_PtVel,
765  const G4FieldTrack& B_PtVel,
766  const G4FieldTrack& SubStart_PtVel,
767  const G4ThreeVector& E_Point,
768  const G4FieldTrack& ApproxIntersecPointV,
769  G4int substep_no, G4int substep_no_p, G4int depth )
770 {
771  // Expect that 'msg' can hold the name of the calling method
772 
773  // FieldTrack 'points' A and B have been tangled
774  // Whereas A should be before B, it is found that curveLen(B) < curveLen(A)
775  G4int verboseLevel= 5;
776  G4double curveDist = B_PtVel.GetCurveLength() - A_PtVel.GetCurveLength();
777  G4VIntersectionLocator::printStatus( A_PtVel, B_PtVel,
778  -1.0, NewSafety, substep_no, msg, verboseLevel );
779  msg << "Error in advancing propagation." << G4endl
780  << " Point A (start) is " << A_PtVel << G4endl
781  << " Point B (end) is " << B_PtVel << G4endl
782  << " Curve distance is " << curveDist << G4endl
783  << G4endl
784  << "The final curve point is not further along"
785  << " than the original!" << G4endl;
786  msg << " Value of fEpsStep= " << epsStep << G4endl;
787 
788  G4int oldprc = msg.precision(20);
789  msg << " Point A (Curve start) is " << StartPointVel << G4endl
790  << " Point B (Curve end) is " << EndPointVel << G4endl
791  << " Point A (Current start) is " << A_PtVel << G4endl
792  << " Point B (Current end) is " << B_PtVel << G4endl
793  << " Point S (Sub start) is " << SubStart_PtVel
794  << " Point E (Trial Point) is " << E_Point << G4endl
795  << " Point F (Intersection) is " << ApproxIntersecPointV
796  << G4endl
797  << " LocateIntersection parameters are : " << G4endl
798  << " Substep no (total) = " << substep_no << G4endl
799  << " Substep (depth= " << depth << substep_no_p;
800  msg.precision(oldprc);
801 }
802 
804 //
805 // ReportProgress.
806 //
807 void G4VIntersectionLocator::ReportProgress( std::ostream& oss,
808  const G4FieldTrack& StartPointVel,
809  const G4FieldTrack& EndPointVel,
810  G4int substep_no,
811  const G4FieldTrack& A_PtVel,
812  const G4FieldTrack& B_PtVel,
813  G4double safetyLast,
814  G4int depth )
815 
816 {
817  oss << "ReportProgress: Current status of intersection search: " << G4endl;
818  if( depth > 0 ) oss << " Depth= " << depth;
819  oss << " Substep no = " << substep_no << G4endl;
820  G4int verboseLevel = 5;
821  G4double safetyPrev = -1.0; // Add as argument ?
822 
823  printStatus( StartPointVel, EndPointVel, -1.0, -1.0, -1,
824  oss, verboseLevel);
825  oss << " * Start and end-point of requested Step:" << G4endl;
826  oss << " ** State of point A: ";
827  printStatus( A_PtVel, A_PtVel, -1.0, safetyPrev, substep_no-1,
828  oss, verboseLevel);
829  oss << " ** State of point B: ";
830  printStatus( A_PtVel, B_PtVel, -1.0, safetyLast, substep_no,
831  oss, verboseLevel);
832 }
833 
835 //
836 // ReportImmediateHit.
837 //
838 void
839 G4VIntersectionLocator::ReportImmediateHit( const char* MethodName,
840  const G4ThreeVector& StartPosition,
841  const G4ThreeVector& TrialPoint,
842  G4double tolerance,
843  unsigned long int numCalls )
844 {
845  static G4ThreadLocal unsigned int occurredOnTop= 0;
846  static G4ThreadLocal G4ThreeVector* ptrLast = nullptr;
847  if( ptrLast == nullptr )
848  {
849  ptrLast= new G4ThreeVector( DBL_MAX, DBL_MAX, DBL_MAX );
850  G4AutoDelete::Register(ptrLast);
851  }
852  G4ThreeVector &lastStart= *ptrLast;
853 
854  if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
855  {
856  static G4ThreadLocal unsigned int numUnmoved = 0;
857  static G4ThreadLocal unsigned int numStill = 0; // Still at same point
858 
859  G4cout << "Intersection F == start A in " << MethodName;
860  G4cout << "Start Point: " << StartPosition << G4endl;
861  G4cout << " Start-Trial: " << TrialPoint - StartPosition;
862  G4cout << " Start-last: " << StartPosition - lastStart;
863 
864  if( (StartPosition - lastStart).mag() < tolerance )
865  {
866  // We are at position of last 'Start' position - ie unmoved
867  ++numUnmoved;
868  ++numStill;
869  G4cout << " { Unmoved: " << " still#= " << numStill
870  << " total # = " << numUnmoved << " } - ";
871  }
872  else
873  {
874  numStill = 0;
875  }
876  G4cout << " Occured: " << ++occurredOnTop;
877  G4cout << " out of total calls= " << numCalls;
878  G4cout << G4endl;
879  lastStart = StartPosition;
880  }
881 } // End of ReportImmediateHit()