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