Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4VIntersectionLocator.cc 110772 2018-06-13 07:47:49Z 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 );
198  auto integrDriver = GetChordFinderFor()->GetIntegrationDriver();
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 //
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  G4Exception("G4VIntersectionLocator::AdjustmentOfFoundIntersection()",
470  "GeomNav0003", JustWarning,
471  "No intersection. Parallels lines!");
472  }
473 #endif
474  lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
475 
476  // New candidate for Intersection
477  //
478  NewPoint = CurrentF_Point+lambda*MomentumDir;
479 
480  // Distance from CurrentF to Calculated Intersection
481  //
482  dist = std::abs(lambda);
483 
484  if ( dist<kCarTolerance*0.001 ) { return false; }
485 
486  // Calculation of new intersection point on the path.
487  //
488  if ( IntersectAF ) // First part intersects
489  {
490  G4double stepLengthFP;
491  G4ThreeVector Point_P = CurrentA_Point;
493  Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
494  fPreviousSafety, fPreviousSftOrigin,
495  stepLengthFP, Point_G );
496 
497  }
498  else // Second part intersects
499  {
500  G4double stepLengthFP;
502  Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
503  fPreviousSafety, fPreviousSftOrigin,
504  stepLengthFP, Point_G );
505  }
506  if ( Intersects_FP )
507  {
508  goodAdjust = true;
509  IntersectionPoint = Point_G;
510  }
511  }
512 
513  return goodAdjust;
514 }
515 
517 //
518 // GetSurfaceNormal.
519 //
521 GetSurfaceNormal(const G4ThreeVector& CurrentInt_Point,
522  G4bool& validNormal) // const
523 {
524  G4ThreeVector NormalAtEntry; // ( -10. , -10., -10. );
525 
526  G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
527  G4bool validNormalLast;
528 
529  // Relies on a call to Navigator::ComputeStep in IntersectChord before
530  // this call
531  //
532  NormalAtEntryLast = GetLastSurfaceNormal( CurrentInt_Point, validNormalLast );
533  // May return valid=false in cases, including
534  // - if the candidate volume was not found (eg exiting world), or
535  // - a replica was involved -- determined the step size.
536  // (This list is not complete.)
537 
538 #ifdef G4DEBUG_FIELD
539  if ( validNormalLast
540  && ( std::fabs(NormalAtEntryLast.mag2() - 1.0) > perThousand ) )
541  {
542  std::ostringstream message;
543  message << "PROBLEM: Normal is not unit - magnitude = "
544  << NormalAtEntryLast.mag() << G4endl;
545  message << " at trial intersection point " << CurrentInt_Point << G4endl;
546  message << " Obtained from Get *Last* Surface Normal.";
547  G4Exception("G4VIntersectionLocator::GetSurfaceNormal()",
548  "GeomNav1002", JustWarning, message);
549  }
550 #endif
551 
552  if( validNormalLast )
553  {
554  NormalAtEntry=NormalAtEntryLast;
555  }
556  validNormal = validNormalLast;
557 
558  return NormalAtEntry;
559 }
560 
562 //
563 // GetGlobalSurfaceNormal.
564 //
566 GetGlobalSurfaceNormal(const G4ThreeVector& CurrentE_Point,
567  G4bool& validNormal)
568 {
569  G4ThreeVector localNormal=
570  GetLocalSurfaceNormal( CurrentE_Point, validNormal );
571  G4AffineTransform localToGlobal= // Must use the same Navigator !!
573  G4ThreeVector globalNormal =
574  localToGlobal.TransformAxis( localNormal );
575 
576 #ifdef G4DEBUG_FIELD
577  if( validNormal && ( std::fabs(globalNormal.mag2() - 1.0) > perThousand ) )
578  {
579  std::ostringstream message;
580  message << "**************************************************************"
581  << G4endl;
582  message << " Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
583  << G4endl;
584  message << " * Constituents: " << G4endl;
585  message << " Local Normal= " << localNormal << G4endl;
586  message << " Transform: " << G4endl
587  << " Net Translation= " << localToGlobal.NetTranslation()
588  << G4endl
589  << " Net Rotation = " << localToGlobal.NetRotation()
590  << G4endl;
591  message << " * Result: " << G4endl;
592  message << " Global Normal= " << localNormal << G4endl;
593  message << "**************************************************************";
594  G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
595  "GeomNav1002", JustWarning, message);
596  }
597 #endif
598 
599  return globalNormal;
600 }
601 
603 //
604 // GetLastSurfaceNormal.
605 //
607 GetLastSurfaceNormal( const G4ThreeVector& intersectPoint,
608  G4bool& normalIsValid) const
609 {
610  G4ThreeVector normalVec;
611  G4bool validNorm;
612  normalVec = fiNavigator->GetGlobalExitNormal( intersectPoint, &validNorm );
613  normalIsValid= validNorm;
614 
615  return normalVec;
616 }
617 
619 //
620 // ReportTrialStep.
621 //
623  const G4ThreeVector& ChordAB_v,
624  const G4ThreeVector& ChordEF_v,
625  const G4ThreeVector& NewMomentumDir,
626  const G4ThreeVector& NormalAtEntry,
627  G4bool validNormal )
628 {
629  G4double ABchord_length = ChordAB_v.mag();
630  G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ) ;
631  G4double MomDir_dot_ABchord;
632  MomDir_dot_ABchord= (1.0 / ABchord_length) * NewMomentumDir.dot( ChordAB_v );
633 
634  std::ostringstream outStream;
635  outStream // G4cout
636  << std::setw(6) << " Step# "
637  << std::setw(17) << " |ChordEF|(mag)" << " "
638  << std::setw(18) << " uMomentum.Normal" << " "
639  << std::setw(18) << " uMomentum.ABdir " << " "
640  << std::setw(16) << " AB-dist " << " "
641  << " Chord Vector (EF) "
642  << G4endl;
643  outStream.precision(7);
644  outStream // G4cout
645  << " " << std::setw(5) << step_no
646  << " " << std::setw(18) << ChordEF_v.mag()
647  << " " << std::setw(18) << MomDir_dot_Norm
648  << " " << std::setw(18) << MomDir_dot_ABchord
649  << " " << std::setw(12) << ABchord_length
650  << " " << ChordEF_v
651  << G4endl;
652  outStream // G4cout
653  << " MomentumDir= " << " " << NewMomentumDir
654  << " Normal at Entry E= " << NormalAtEntry
655  << " AB chord = " << ChordAB_v
656  << G4endl;
657  G4cout << outStream.str(); // ostr_verbose;
658 
659  if( ( std::fabs(NormalAtEntry.mag2() - 1.0) > perThousand ) )
660  {
661  std::ostringstream message;
662  message << "Normal is not unit - mag= " << NormalAtEntry.mag() << G4endl
663  << " ValidNormalAtE = " << validNormal;
664  G4Exception("G4VIntersectionLocator::ReportTrialStep()",
665  "GeomNav1002", JustWarning, message);
666  }
667  return;
668 }
669 
671 //
672 // LocateGlobalPointWithinVolumeAndCheck.
673 //
674 // Locate point using navigator: updates state of Navigator
675 // By default, it assumes that the point is inside the current volume,
676 // and returns true.
677 // In check mode, checks whether the point is *inside* the volume.
678 // If it is inside, it returns true
679 // If not, issues a warning and returns false.
680 //
683 {
684  G4bool good= true;
686  const G4String
687  MethodName("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()");
688 
689  if( fCheckMode )
690  {
691  G4bool navCheck= nav->IsCheckModeActive(); // Recover original value
692  nav->CheckMode(true);
693 
694  // Identify the current volume
695 
697  G4VPhysicalVolume* motherPhys= startTH->GetVolume();
698  G4VSolid* motherSolid= startTH->GetSolid();
699  G4AffineTransform transform = nav->GetGlobalToLocalTransform();
700  // GetLocalToGlobalTransform();
701  G4int motherCopyNo= motherPhys->GetCopyNo();
702 
703  // Let's check that the point is inside the current solid
704  G4ThreeVector localPosition = transform.TransformPoint(position);
705  EInside inMother= motherSolid->Inside( localPosition );
706  if( inMother != kInside )
707  {
708  std::ostringstream message;
709  message << "Position located "
710  << ( inMother == kSurface ? " on Surface " : " outside " )
711  << "expected volume" << G4endl
712  << " Safety (from Outside) = "
713  << motherSolid->DistanceToIn(localPosition);
714  G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
715  "GeomNav1002", JustWarning, message);
716  }
717 
718  // 1. Simple next step - quick relocation and check result.
719  // nav->LocateGlobalPointWithinVolume( position );
720 
721  // 2. Full relocation - to cross-check answer !
722  G4VPhysicalVolume* nextPhysical= nav->LocateGlobalPointAndSetup(position);
723  if( (nextPhysical != motherPhys)
724  || (nextPhysical->GetCopyNo() != motherCopyNo )
725  )
726  {
727  G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
728  "GeomNav1002", JustWarning,
729  "Position located outside expected volume.");
730  }
731  nav->CheckMode(navCheck); // Recover original value
732  }
733  else
734  {
735  nav->LocateGlobalPointWithinVolume( position );
736  }
737  return good;
738 }
739 
741 //
742 // LocateGlobalPointWithinVolumeCheckAndReport.
743 //
746  const G4String& CodeLocationInfo,
747  G4int /* CheckMode */)
748 {
749  // Save value of Check mode first
750  G4bool oldCheck= GetCheckMode();
751 
753  if( !ok )
754  {
755  std::ostringstream message;
756  message << "Failed point location." << G4endl
757  << " Code Location info: " << CodeLocationInfo;
758  G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport()",
759  "GeomNav1002", JustWarning, message);
760 /*
761  if( CheckMode > 1 )
762  {
763  // Additional information
764  }
765 */
766  }
767 
768  SetCheckMode( oldCheck );
769 }
770 
772 //
773 // ReportReversedPoints.
774 //
776 ReportReversedPoints( std::ostringstream& msg,
777  const G4FieldTrack& StartPointVel,
778  const G4FieldTrack& EndPointVel,
779  G4double NewSafety, G4double epsStep,
780  const G4FieldTrack& A_PtVel,
781  const G4FieldTrack& B_PtVel,
782  const G4FieldTrack& SubStart_PtVel,
783  const G4ThreeVector& E_Point,
784  const G4FieldTrack& ApproxIntersecPointV,
785  G4int substep_no, G4int substep_no_p, G4int depth )
786 {
787  // Expect that 'msg' can hold the name of the calling method
788 
789  // FieldTrack 'points' A and B have been tangled
790  // Whereas A should be before B, it is found that curveLen(B) < curveLen(A)
791  G4int verboseLevel= 5;
792  G4double curveDist = B_PtVel.GetCurveLength() - A_PtVel.GetCurveLength();
793  G4VIntersectionLocator::printStatus( A_PtVel, B_PtVel,
794  -1.0, NewSafety, substep_no, msg, verboseLevel );
795  msg << "Error in advancing propagation." << G4endl
796  << " Point A (start) is " << A_PtVel << G4endl
797  << " Point B (end) is " << B_PtVel << G4endl
798  << " Curve distance is " << curveDist << G4endl
799  << G4endl
800  << "The final curve point is not further along"
801  << " than the original!" << G4endl;
802  msg << " Value of fEpsStep= " << epsStep << G4endl;
803 
804  G4int oldprc = msg.precision(20);
805  msg << " Point A (Curve start) is " << StartPointVel << G4endl
806  << " Point B (Curve end) is " << EndPointVel << G4endl
807  << " Point A (Current start) is " << A_PtVel << G4endl
808  << " Point B (Current end) is " << B_PtVel << G4endl
809  << " Point S (Sub start) is " << SubStart_PtVel
810  << " Point E (Trial Point) is " << E_Point << G4endl
811  << " Point F (Intersection) is " << ApproxIntersecPointV
812  << G4endl
813  << " LocateIntersection parameters are : " << G4endl
814  << " Substep no (total) = " << substep_no << G4endl
815  << " Substep (depth= " << depth << substep_no_p;
816  msg.precision(oldprc);
817 }
818 
820 //
821 // ReportProgress.
822 //
824  const G4FieldTrack& StartPointVel,
825  const G4FieldTrack& EndPointVel,
826  G4int substep_no,
827  const G4FieldTrack& A_PtVel,
828  const G4FieldTrack& B_PtVel,
829  G4double safetyLast,
830  G4int depth )
831 
832 {
833  oss << "ReportProgress: Current status of intersection search: " << G4endl;
834  if( depth > 0 ) oss << " Depth= " << depth;
835  oss << " Substep no = " << substep_no << G4endl;
836  G4int verboseLevel = 5;
837 
838  // printStatus args: (FT0, FT1, dRequestStep, dSafety, iStepNum, os, iVerb);
839  G4double safetyPrev = -1.0; // Add as argument ?
840 
841  printStatus( StartPointVel, EndPointVel, -1.0, -1.0, -1,
842  oss, verboseLevel);
843  oss << " * Start and end-point of requested Step:" << G4endl;
844  oss << " ** State of point A: ";
845  printStatus( A_PtVel, A_PtVel, -1.0, safetyPrev, substep_no-1,
846  oss, verboseLevel);
847  oss << " ** State of point B: ";
848  printStatus( A_PtVel, B_PtVel, -1.0, safetyLast, substep_no,
849  oss, verboseLevel);
850 }
851 
853 //
854 // ReportImmediateHit.
855 //
856 void
858  const G4ThreeVector& StartPosition,
859  const G4ThreeVector& TrialPoint,
860  G4double tolerance,
861  unsigned long int numCalls )
862 {
863  static G4ThreadLocal unsigned int occurredOnTop= 0;
864  static G4ThreadLocal G4ThreeVector *ptrLast= 0;
865  if( !ptrLast )
866  {
867  ptrLast= new G4ThreeVector( DBL_MAX, DBL_MAX, DBL_MAX );
868  G4AutoDelete::Register(ptrLast);
869  }
870  G4ThreeVector &lastStart= *ptrLast;
871 
872  if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
873  {
874  static G4ThreadLocal unsigned int numUnmoved= 0;
875  static G4ThreadLocal unsigned int numStill= 0; // Still at same point
876 
877  G4cout << "Intersection F == start A in " << MethodName;
878  G4cout << "Start Point: " << StartPosition << G4endl;
879  // G4cout << "Trial Point: " << TrialPoint << G4endl;
880  G4cout << " Start-Trial: " << TrialPoint - StartPosition; // << G4endl;
881  G4cout << " Start-last: " << StartPosition - lastStart;
882 
883  if( (StartPosition - lastStart).mag() < tolerance )
884  {
885  // We are at position of last 'Start' position - ie unmoved
886  ++numUnmoved;
887  ++numStill;
888  G4cout << " { Unmoved: " << " still#= " << numStill
889  << " total # = " << numUnmoved << " } - ";
890  }
891  else
892  {
893  numStill = 0;
894  }
895  G4cout << " Occured: " << ++occurredOnTop;
896  G4cout << " out of total calls= " << numCalls;
897  G4cout << G4endl;
898  lastStart = StartPosition;
899  }
900 } // End of ReportImmediateHit()
G4VIntegrationDriver * GetIntegrationDriver()
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
static constexpr double perThousand
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetLocalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
void LocateGlobalPointWithinVolumeCheckAndReport(const G4ThreeVector &pos, const G4String &CodeLocationInfo, G4int CheckMode)
G4LogicalVolume * GetLogicalVolume() const
G4bool IsCheckModeActive() const
G4bool LocateGlobalPointWithinVolumeAndCheck(const G4ThreeVector &pos)
G4ChordFinder * GetChordFinderFor()
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:588
#define G4endl
Definition: G4ios.hh:61
G4TouchableHistory * CreateTouchableHistory() const
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
#define fPreviousSafety
double z() const
double dot(const Hep3Vector &) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
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)
#define G4ThreadLocal
Definition: tls.hh:69
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4ThreeVector GetLastSurfaceNormal(const G4ThreeVector &intersectPoint, G4bool &validNormal) const
void Register(T *inst)
Definition: G4AutoDelete.hh:65
void CheckMode(G4bool mode)
virtual EInside Inside(const G4ThreeVector &p) const =0
const G4AffineTransform & GetGlobalToLocalTransform() const
const G4AffineTransform GetLocalToGlobalTransform() const
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4ThreeVector GetPosition() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4VSolid * GetSolid(G4int depth=0) const
G4ThreeVector GetGlobalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
const G4NavigationHistory * GetHistory() const
#define fPreviousSftOrigin
G4VPhysicalVolume * GetVolume(G4int depth=0) const
G4Navigator * GetNavigatorFor()
void SetWorldVolume(G4VPhysicalVolume *pWorld)
const G4AffineTransform & GetTopTransform() const
G4ThreeVector GetMomentum() const
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
double mag2() const
virtual G4int GetCopyNo() const =0
G4GLOB_DLL std::ostream G4cerr
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)
G4ThreeVector NetTranslation() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
double mag() const
virtual G4TouchableHistoryHandle CreateTouchableHistoryHandle() const
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
G4TouchableHistory * fpTouchable
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
const G4ThreeVector & GetMomentumDir() const
G4GLOB_DLL std::ostream G4cout
double x() const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4VSolid * GetSolid() const
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:131
static G4GeometryTolerance * GetInstance()
double y() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetEpsilonStepFor()
#define DBL_MAX
Definition: templates.hh:83
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)
static constexpr double perThousand
Definition: G4SIunits.hh:333
G4double GetSurfaceTolerance() const
G4double GetCurveLength() const
G4VIntersectionLocator(G4Navigator *theNavigator)
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
G4RotationMatrix NetRotation() const
virtual G4ThreeVector GetGlobalExitNormal(const G4ThreeVector &point, G4bool *valid)