Geant4-11
Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | Static Private Attributes
G4MultiLevelLocator Class Reference

#include <G4MultiLevelLocator.hh>

Inheritance diagram for G4MultiLevelLocator:
G4VIntersectionLocator

Public Member Functions

void AddAdjustementOfFoundIntersection (G4bool UseCorrection)
 
void AdjustIntersections (G4bool UseCorrection)
 
G4bool AreIntersectionsAdjusted ()
 
G4bool EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
 
 G4MultiLevelLocator (G4Navigator *theNavigator)
 
G4bool GetAdjustementOfFoundIntersection ()
 
G4bool GetCheckMode ()
 
G4ChordFinderGetChordFinderFor ()
 
G4double GetDeltaIntersectionFor ()
 
G4double GetEpsilonStepFor ()
 
G4NavigatorGetNavigatorFor ()
 
G4int GetVerboseFor ()
 
G4bool IntersectChord (const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=nullptr)
 
void printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
 
void ReportStatistics ()
 
void SetCheckMode (G4bool value)
 
void SetChordFinderFor (G4ChordFinder *fCFinder)
 
void SetDeltaIntersectionFor (G4double deltaIntersection)
 
void SetEpsilonStepFor (G4double EpsilonStep)
 
void SetMaxSteps (unsigned int valMax)
 
void SetNavigatorFor (G4Navigator *fNavigator)
 
void SetSafetyParametersFor (G4bool UseSafety)
 
void SetVerboseFor (G4int fVerbose)
 
void SetWarnSteps (unsigned int valWarn)
 
 ~G4MultiLevelLocator ()
 

Static Public Member Functions

static void printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum, std::ostream &oss, G4int verboseLevel)
 

Protected Member Functions

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)
 
G4bool CheckAndReEstimateEndpoint (const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
 
G4ThreeVector GetGlobalSurfaceNormal (const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
 
G4ThreeVector GetSurfaceNormal (const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
 
G4bool LocateGlobalPointWithinVolumeAndCheck (const G4ThreeVector &pos)
 
void LocateGlobalPointWithinVolumeCheckAndReport (const G4ThreeVector &pos, const G4String &CodeLocationInfo, G4int CheckMode)
 
G4FieldTrack ReEstimateEndpoint (const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
 
void ReportImmediateHit (const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, G4double tolerance, unsigned long int numCalls)
 
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)
 
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)
 

Protected Attributes

G4bool fCheckMode = false
 
G4NavigatorfHelpingNavigator
 
G4ChordFinderfiChordFinder = nullptr
 
G4double fiDeltaIntersection = -1.0
 
G4double fiEpsilonStep = -1.0
 
G4NavigatorfiNavigator
 
G4bool fiUseSafety = false
 
G4TouchableHistoryfpTouchable = nullptr
 
G4bool fUseNormalCorrection = false
 
G4int fVerboseLevel = 0
 
G4double kCarTolerance
 

Private Member Functions

G4ThreeVector GetLastSurfaceNormal (const G4ThreeVector &intersectPoint, G4bool &validNormal) const
 
G4ThreeVector GetLocalSurfaceNormal (const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
 
void ReportFieldValue (const G4FieldTrack &locationPV, const char *nameLoc, const G4EquationOfMotion *equation)
 

Private Attributes

unsigned int fMaxSteps = 10000
 
unsigned long int fNumAdvanceFull = 0
 
unsigned long int fNumAdvanceGood = 0
 
unsigned long int fNumAdvanceTrials = 0
 
unsigned long int fNumCalls = 0
 
unsigned int fWarnSteps = 1000
 
G4FieldTrackptrInterMedFT [max_depth+1]
 

Static Private Attributes

static const G4int max_depth = 10
 

Detailed Description

Definition at line 47 of file G4MultiLevelLocator.hh.

Constructor & Destructor Documentation

◆ G4MultiLevelLocator()

G4MultiLevelLocator::G4MultiLevelLocator ( G4Navigator theNavigator)

Definition at line 39 of file G4MultiLevelLocator.cc.

40 : G4VIntersectionLocator(theNavigator)
41{
42 // In case of too slow progress in finding Intersection Point
43 // intermediates Points on the Track must be stored.
44 // Initialise the array of Pointers [max_depth+1] to do this
45
46 G4ThreeVector zeroV(0.0,0.0,0.0);
47 for ( auto idepth=0; idepth<max_depth+1; ++idepth )
48 {
49 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
50 }
51
52 if (fCheckMode)
53 {
54 // Trial values Loose Medium Tight
55 // To happen: Infrequent Occasional Often
56 SetMaxSteps(150); // 300 85 25
57 SetWarnSteps(80); // 250 60 15
58 }
59}
G4FieldTrack * ptrInterMedFT[max_depth+1]
static const G4int max_depth
void SetMaxSteps(unsigned int valMax)
void SetWarnSteps(unsigned int valWarn)
G4VIntersectionLocator(G4Navigator *theNavigator)

References G4VIntersectionLocator::fCheckMode, max_depth, ptrInterMedFT, SetMaxSteps(), and SetWarnSteps().

◆ ~G4MultiLevelLocator()

G4MultiLevelLocator::~G4MultiLevelLocator ( )

Definition at line 61 of file G4MultiLevelLocator.cc.

62{
63 for ( auto idepth=0; idepth<max_depth+1; ++idepth )
64 {
65 delete ptrInterMedFT[idepth];
66 }
67#ifdef G4DEBUG_FIELD
69#endif
70}

References max_depth, ptrInterMedFT, and ReportStatistics().

Member Function Documentation

◆ AddAdjustementOfFoundIntersection()

void G4VIntersectionLocator::AddAdjustementOfFoundIntersection ( G4bool  UseCorrection)
inlineinherited

◆ AdjustIntersections()

void G4VIntersectionLocator::AdjustIntersections ( G4bool  UseCorrection)
inlineinherited

◆ AdjustmentOfFoundIntersection()

G4bool G4VIntersectionLocator::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 
)
protectedinherited

Definition at line 439 of file G4VIntersectionLocator.cc.

449{
450 G4double dist,lambda;
451 G4ThreeVector Normal, NewPoint, Point_G;
452 G4bool goodAdjust = false, Intersects_FP = false, validNormal = false;
453
454 // Get SurfaceNormal of Intersecting Solid
455 //
456 Normal = GetGlobalSurfaceNormal(CurrentE_Point,validNormal);
457 if(!validNormal) { return false; }
458
459 // Intersection between Line and Plane
460 //
461 G4double n_d_m = Normal.dot(MomentumDir);
462 if ( std::abs(n_d_m)>kCarTolerance )
463 {
464#ifdef G4VERBOSE
465 if ( fVerboseLevel>1 )
466 {
467 G4Exception("G4VIntersectionLocator::AdjustmentOfFoundIntersection()",
468 "GeomNav0003", JustWarning,
469 "No intersection. Parallels lines!");
470 }
471#endif
472 lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
473
474 // New candidate for Intersection
475 //
476 NewPoint = CurrentF_Point+lambda*MomentumDir;
477
478 // Distance from CurrentF to Calculated Intersection
479 //
480 dist = std::abs(lambda);
481
482 if ( dist<kCarTolerance*0.001 ) { return false; }
483
484 // Calculation of new intersection point on the path.
485 //
486 if ( IntersectAF ) // First part intersects
487 {
488 G4double stepLengthFP;
489 G4ThreeVector Point_P = CurrentA_Point;
491 Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
493 stepLengthFP, Point_G );
494
495 }
496 else // Second part intersects
497 {
498 G4double stepLengthFP;
500 Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
502 stepLengthFP, Point_G );
503 }
504 if ( Intersects_FP )
505 {
506 goodAdjust = true;
507 IntersectionPoint = Point_G;
508 }
509 }
510
511 return goodAdjust;
512}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
#define fPreviousSftOrigin
#define fPreviousSafety
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
double dot(const Hep3Vector &) const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:614
G4ThreeVector GetGlobalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
G4Navigator * GetNavigatorFor()
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=nullptr)

References CLHEP::Hep3Vector::dot(), fPreviousSafety, fPreviousSftOrigin, G4VIntersectionLocator::fVerboseLevel, G4Exception(), G4VIntersectionLocator::GetGlobalSurfaceNormal(), G4VIntersectionLocator::GetNavigatorFor(), G4VIntersectionLocator::IntersectChord(), JustWarning, G4VIntersectionLocator::kCarTolerance, G4InuclParticleNames::lambda, and G4Navigator::LocateGlobalPointWithinVolume().

Referenced by G4BrentLocator::EstimateIntersectionPoint(), EstimateIntersectionPoint(), and G4SimpleLocator::EstimateIntersectionPoint().

◆ AreIntersectionsAdjusted()

G4bool G4VIntersectionLocator::AreIntersectionsAdjusted ( )
inlineinherited

◆ CheckAndReEstimateEndpoint()

G4bool G4VIntersectionLocator::CheckAndReEstimateEndpoint ( const G4FieldTrack CurrentStartA,
const G4FieldTrack EstimatedEndB,
G4FieldTrack RevisedEndPoint,
G4int errorCode 
)
protectedinherited

Definition at line 329 of file G4VIntersectionLocator.cc.

334{
335 G4double linDistSq, curveDist;
336
337 G4bool recalculated = false;
338 curveError= 0;
339
340 linDistSq = ( EstimatedEndB.GetPosition()
341 - CurrentStartA.GetPosition() ).mag2();
342 curveDist = EstimatedEndB.GetCurveLength()
343 - CurrentStartA.GetCurveLength();
344 if( (curveDist>=0.0)
345 && (curveDist*curveDist *(1.0+2.0*fiEpsilonStep ) < linDistSq ) )
346 {
347 G4FieldTrack newEndPointFT = EstimatedEndB; // Unused
348
349 if (curveDist>0.0)
350 {
351 // Re-integrate to obtain a new B
352 RevisedEndPoint = ReEstimateEndpoint( CurrentStartA,
353 EstimatedEndB,
354 linDistSq,
355 curveDist );
356 recalculated = true;
357 }
358 else
359 {
360 // Zero length -> no advance!
361 newEndPointFT = CurrentStartA;
362 recalculated = true;
363 curveError = 1; // Unexpected co-incidence - milder mixup
364
365 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
366 "GeomNav1002", JustWarning,
367 "A & B are at equal distance in 2nd half. A & B will coincide." );
368 }
369 }
370
371 // Sanity check
372 //
373 if( curveDist < 0.0 )
374 {
375 curveError = 2; // Real mixup
376 }
377 return recalculated;
378}
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)

References G4VIntersectionLocator::fiEpsilonStep, G4Exception(), G4FieldTrack::GetCurveLength(), G4FieldTrack::GetPosition(), JustWarning, and G4VIntersectionLocator::ReEstimateEndpoint().

Referenced by EstimateIntersectionPoint().

◆ EstimateIntersectionPoint()

G4bool G4MultiLevelLocator::EstimateIntersectionPoint ( const G4FieldTrack curveStartPointTangent,
const G4FieldTrack curveEndPointTangent,
const G4ThreeVector trialPoint,
G4FieldTrack intersectPointTangent,
G4bool recalculatedEndPoint,
G4double fPreviousSafety,
G4ThreeVector fPreviousSftOrigin 
)
virtual

Implements G4VIntersectionLocator.

Definition at line 115 of file G4MultiLevelLocator.cc.

123{
124 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
125 const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()";
126
127 G4bool found_approximate_intersection = false;
128 G4bool there_is_no_intersection = false;
129
130 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
131 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
132 G4ThreeVector CurrentE_Point = TrialPoint;
133 G4bool validNormalAtE = false;
134 G4ThreeVector NormalAtEntry;
135
136 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
137 G4bool validIntersectP= true; // Is it current ?
138 G4double NewSafety = 0.0;
139 G4bool last_AF_intersection = false;
140
141 auto integrDriver = GetChordFinderFor()->GetIntegrationDriver();
142 G4bool driverReIntegrates = integrDriver->DoesReIntegrate();
143
144 G4bool first_section = true;
145 recalculatedEndPoint = false;
146 G4bool restoredFullEndpoint = false;
147
148 unsigned int substep_no = 0;
149
150 // Statistics for substeps
151 static G4ThreadLocal unsigned int max_no_seen= 0;
152
153#ifdef G4DEBUG_FIELD
154 unsigned int trigger_substepno_print = 0;
155 const G4double tolerance = 1.0e-8 * CLHEP::mm;
156 unsigned int biggest_depth = 0;
157 // using kInitialisingCL = G4LocatorChangeRecord::kInitialisingCL;
158#endif
159
160 // Log the location, iteration of changes in A,B
161 //----------------------------------------------
162 static G4ThreadLocal G4LocatorChangeLogger endChangeA("StartPointA"),
163 endChangeB("EndPointB"), recApproxPoint("ApproxPoint"),
164 pointH_logger("Trial points 'E': position, normal");
165 unsigned int eventCount = 0;
166
167 if (fCheckMode)
168 {
169 // Clear previous call's data
170 endChangeA.clear();
171 endChangeB.clear();
172 recApproxPoint.clear();
173 pointH_logger.clear();
174
175 // Record the initialisation
176 ++eventCount;
177 endChangeA.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no,
178 eventCount, CurrentA_PointVelocity );
179 endChangeB.AddRecord( G4LocatorChangeRecord::kInitialisingCL, substep_no,
180 eventCount, CurrentB_PointVelocity );
181 }
182
183 //--------------------------------------------------------------------------
184 // Algorithm for the case if progress in founding intersection is too slow.
185 // Process is defined too slow if after N=param_substeps advances on the
186 // path, it will be only 'fraction_done' of the total length.
187 // In this case the remaining length is divided in two half and
188 // the loop is restarted for each half.
189 // If progress is still too slow, the division in two halfs continue
190 // until 'max_depth'.
191 //--------------------------------------------------------------------------
192
193 const G4int param_substeps = 5; // Test value for the maximum number
194 // of substeps
195 const G4double fraction_done = 0.3;
196
197 G4bool Second_half = false; // First half or second half of divided step
198
199 // We need to know this for the 'final_section':
200 // real 'final_section' or first half 'final_section'
201 // In algorithm it is considered that the 'Second_half' is true
202 // and it becomes false only if we are in the first-half of level
203 // depthness or if we are in the first section
204
205 unsigned int depth = 0; // Depth counts subdivisions of initial step made
206 ++fNumCalls;
207
208 NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
209
210 if (fCheckMode)
211 {
212 pointH_logger.AddRecord( G4LocatorChangeRecord::kInitialisingCL,
213 substep_no, eventCount,
214 G4FieldTrack( CurrentE_Point,0.,NormalAtEntry,0.,
215 0., 1.,G4ThreeVector(0.),0.,0.) );
216 #if (G4DEBUG_FIELD>1)
217 G4ThreeVector StartPosition = CurveStartPointVelocity.GetPosition();
218 if( (TrialPoint - StartPosition).mag2() < sqr(tolerance))
219 {
220 ReportImmediateHit( MethodName, StartPosition, TrialPoint,
221 tolerance, fNumCalls);
222 }
223 #endif
224 }
225
226 // Intermediates Points on the Track = Subdivided Points must be stored.
227 // Give the initial values to 'InterMedFt'
228 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
229 //
230 *ptrInterMedFT[0] = CurveEndPointVelocity;
231 for ( auto idepth=1; idepth<max_depth+1; ++idepth )
232 {
233 *ptrInterMedFT[idepth] = CurveStartPointVelocity;
234 }
235
236 // Final_section boolean store
237 //
238 G4bool fin_section_depth[max_depth];
239 for ( auto idepth=0; idepth<max_depth; ++idepth )
240 {
241 fin_section_depth[idepth] = true;
242 }
243 // 'SubStartPoint' is needed to calculate the length of the divided step
244 //
245 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
246
247 do // Loop checking, 07.10.2016, J.Apostolakis
248 {
249 unsigned int substep_no_p = 0;
250 G4bool sub_final_section = false; // the same as final_section,
251 // but for 'sub_section'
252 SubStart_PointVelocity = CurrentA_PointVelocity;
253
254 do // Loop checking, 07.10.2016, J.Apostolakis
255 { // REPEAT param
256
257#ifdef G4DEBUG_FIELD
258 if( CurrentA_PointVelocity.GetCurveLength() >=
259 CurrentB_PointVelocity.GetCurveLength() )
260 {
261 G4cerr << "ERROR> (Start) Point A coincides with or has gone past (end) point B"
262 << "MLL: iters = " << substep_no << G4endl;
263 // G4LocatorChangeRecord::ReportVector(G4cerr, "endPointB", endChangeB );
264 // G4cerr<<"EndPoints A(start) and B(end): combined changes " << G4endl;
265 G4LocatorChangeLogger::ReportEndChanges(G4cerr, endChangeA, endChangeB);
266 }
267#endif
268 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
269 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
270
271 // F = a point on true AB path close to point E
272 // (the closest if possible)
273 //
274 ApproxIntersecPointV = GetChordFinderFor()
275 ->ApproxCurvePointV( CurrentA_PointVelocity,
276 CurrentB_PointVelocity,
277 CurrentE_Point,
279 // The above method is the key & most intuitive part ...
280
281#ifdef G4DEBUG_FIELD
283 substep_no, eventCount, ApproxIntersecPointV ) );
284 G4double lenIntsc= ApproxIntersecPointV.GetCurveLength();
285 G4double lenB = CurrentB_PointVelocity.GetCurveLength();
286 G4double checkVsEnd= lenB - lenIntsc;
287
288 if( lenIntsc > lenB )
289 {
290 std::ostringstream errmsg;
291 errmsg.precision(17);
292 G4double ratio = checkVsEnd / lenB;
293 G4double ratioTol = std::fabs(ratio) / tolerance;
294 errmsg << "Intermediate F point is past end B point" << G4endl
295 << " l( intersection ) = " << lenIntsc << G4endl
296 << " l( endpoint ) = " << lenB << G4endl;
297 errmsg.precision(8);
298 errmsg << " l_end - l_inters = " << checkVsEnd << G4endl
299 << " / l_end = " << ratio << G4endl
300 << " ratio / tolerance = " << ratioTol << G4endl;
301 if( ratioTol < 1.0 )
302 G4Exception(MethodName, "GeomNav0003", JustWarning, errmsg );
303 else
304 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg );
305 }
306#endif
307
308 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
309
310 // First check whether EF is small - then F is a good approx. point
311 // Calculate the length and direction of the chord AF
312 //
313 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
314
315 G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
316 G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry );
317
318#ifdef G4DEBUG_FIELD
319 if( fVerboseLevel > 3 )
320 {
321 G4ThreeVector ChordAB = Point_B - Point_A;
322 G4double ABchord_length = ChordAB.mag();
323 G4double MomDir_dot_ABchord;
324 MomDir_dot_ABchord = (1.0 / ABchord_length)
325 * NewMomentumDir.dot( ChordAB );
326 G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
327 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
328 G4cout << " dot( MomentumDir, ABchord_unit ) = "
329 << MomDir_dot_ABchord << G4endl;
330 }
331#endif
332 G4bool adequate_angle =
333 ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
334 || (! validNormalAtE ); // Invalid, cannot use this criterion
335 G4double EF_dist2 = ChordEF_Vector.mag2();
336 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
337 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
338 {
339 found_approximate_intersection = true;
340
341 // Create the "point" return value
342 //
343 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
344 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
345
347 {
348 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
349 //
350 G4ThreeVector IP;
351 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
352 G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
353 CurrentE_Point, CurrentF_Point, MomentumDir,
354 last_AF_intersection, IP, NewSafety,
355 previousSafety, previousSftOrigin );
356 if ( goodCorrection )
357 {
358 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
359 IntersectedOrRecalculatedFT.SetPosition(IP);
360 }
361 }
362 // Note: in order to return a point on the boundary,
363 // we must return E. But it is F on the curve.
364 // So we must "cheat": we are using the position at point E
365 // and the velocity at point F !!!
366 //
367 // This must limit the length we can allow for displacement!
368 }
369 else // E is NOT close enough to the curve (ie point F)
370 {
371 // Check whether any volumes are encountered by the chord AF
372 // ---------------------------------------------------------
373 // First relocate to restore any Voxel etc information
374 // in the Navigator before calling ComputeStep()
375 //
377
378 G4ThreeVector PointG; // Candidate intersection point
379 G4double stepLengthAF;
380 G4bool Intersects_FB = false;
381 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
382 NewSafety, previousSafety,
383 previousSftOrigin,
384 stepLengthAF,
385 PointG );
386 last_AF_intersection = Intersects_AF;
387 if( Intersects_AF )
388 {
389 // G is our new Candidate for the intersection point.
390 // It replaces "E" and we will see if it's good enough.
391 CurrentB_PointVelocity = ApproxIntersecPointV; // B <- F
392 CurrentE_Point = PointG; // E <- G
393
394 validIntersectP = true; // 'E' has been updated.
395
396 G4bool validNormalLast;
397 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
398 validNormalAtE = validNormalLast;
399
400 // As we move point B, must take care in case the current
401 // AF has no intersection to try current FB!!
402 fin_section_depth[depth] = false;
403
404 if (fCheckMode)
405 {
406 ++eventCount;
407 endChangeB.push_back(
409 substep_no, eventCount, CurrentB_PointVelocity) );
410 }
411#ifdef G4VERBOSE
412 if( fVerboseLevel > 3 )
413 {
414 G4cout << "G4PiF::LI> Investigating intermediate point"
415 << " at s=" << ApproxIntersecPointV.GetCurveLength()
416 << " on way to full s="
417 << CurveEndPointVelocity.GetCurveLength() << G4endl;
418 }
419#endif
420 }
421 else // not Intersects_AF
422 {
423 // In this case:
424 // There is NO intersection of AF with a volume boundary.
425 // We must continue the search in the segment FB!
426 //
428
429 G4double stepLengthFB;
430 G4ThreeVector PointH;
431
432 // Check whether any volumes are encountered by the chord FB
433 // ---------------------------------------------------------
434
435 Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
436 NewSafety, previousSafety,
437 previousSftOrigin,
438 stepLengthFB,
439 PointH );
440 if( Intersects_FB )
441 {
442 // There is an intersection of FB with a volume boundary
443 // H <- First Intersection of Chord FB
444
445 // H is our new Candidate for the intersection point.
446 // It replaces "E" and we will repeat the test to see if
447 // it is a good enough approximate point for us.
448
449 // Note that F must be in volume volA (the same as A)
450 // (otherwise AF would meet a volume boundary!)
451 // A <- F
452 // E <- H
453 //
454 CurrentA_PointVelocity = ApproxIntersecPointV;
455 CurrentE_Point = PointH;
456
457 validIntersectP = true; // 'E' has been updated.
458
459 G4bool validNormalH;
460 NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
461 validNormalAtE = validNormalH;
462
463 if (fCheckMode)
464 {
465 ++eventCount;
466 endChangeA.push_back(
468 substep_no, eventCount, CurrentA_PointVelocity) );
469 G4FieldTrack intersectH_pn('0'); // Point and normal
470 // nothing else will be valid
471 intersectH_pn.SetPosition( PointH );
472 intersectH_pn.SetMomentum( NormalAtEntry );
473 pointH_logger.AddRecord(G4LocatorChangeRecord::kIntersectsFB,
474 substep_no, eventCount, intersectH_pn );
475 }
476 }
477 else // not Intersects_FB
478 {
479 if( fin_section_depth[depth] )
480 {
481 // If B is the original endpoint, this means that whatever
482 // volume(s) intersected the original chord, none touch the
483 // smaller chords we have used.
484 // The value of 'IntersectedOrRecalculatedFT' returned is
485 // likely not valid
486
487 // Check on real final_section or SubEndSection
488 //
489 if( ((Second_half)&&(depth==0)) || (first_section) )
490 {
491 there_is_no_intersection = true; // real final_section
492 }
493 else
494 {
495 // end of subsection, not real final section
496 // exit from the and go to the depth-1 level
497 substep_no_p = param_substeps+2; // exit from the loop
498
499 // but 'Second_half' is still true because we need to find
500 // the 'CurrentE_point' for the next loop
501 Second_half = true;
502 sub_final_section = true;
503 }
504 }
505 else
506 {
507 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
508 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
509 : *ptrInterMedFT[depth] ;
510 SubStart_PointVelocity = CurrentA_PointVelocity;
511 restoredFullEndpoint = true;
512
513 validIntersectP = false; // 'E' has NOT been updated.
514
515 if (fCheckMode)
516 {
517 ++eventCount;
518 endChangeA.push_back(
521 substep_no, eventCount, CurrentA_PointVelocity) );
522 endChangeB.push_back(
525 substep_no, eventCount, CurrentB_PointVelocity) );
526 }
527 }
528 } // Endif (Intersects_FB)
529 } // Endif (Intersects_AF)
530
531 G4int errorEndPt = 0; // Default: no error (if not calling CheckAnd...
532
533 G4bool recalculatedB= false;
534 if( driverReIntegrates )
535 {
536 G4FieldTrack RevisedB_FT = CurrentB_PointVelocity;
537 recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity,
538 CurrentB_PointVelocity,
539 RevisedB_FT,
540 errorEndPt );
541 if( recalculatedB )
542 {
543 CurrentB_PointVelocity = RevisedB_FT; // Use it !
544 // Do not invalidate intersection F -- it is still roughly OK.
545 //
546 // The best course would be to invalidate (reset validIntersectP)
547 // BUT if we invalidate it, we must re-estimate it somewhere! E.g.
548 // validIntersectP = false; // 'E' has NOT been updated.
549
550 if ( (fin_section_depth[depth]) // real final section
551 &&( first_section || ((Second_half)&&(depth==0)) ) )
552 {
553 recalculatedEndPoint = true;
554 IntersectedOrRecalculatedFT = RevisedB_FT;
555 // So that we can return it, if it is the endpoint!
556 }
557 // else
558 // Move forward the other points
559 // - or better flag it, so that they are re-computed when next used
560 // [ Implementation: a counter for # of recomputations
561 // => avoids extra work]
562 }
563 // else
564 // Move forward the other points
565 // - or better flag it, so that they are re-computed when next used
566 // [ Implementation: a counter for # of recomputations
567 // => avoids extra work]
568 if (fCheckMode)
569 {
570 ++eventCount;
571 endChangeB.push_back(
573 substep_no, eventCount, RevisedB_FT ) );
574 }
575 }
576 else
577 {
578 if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() )
579 errorEndPt = 2;
580 }
581
582 if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
583 {
584 std::ostringstream errmsg;
585 ReportReversedPoints(errmsg,
586 CurveStartPointVelocity, CurveEndPointVelocity,
587 NewSafety, fiEpsilonStep,
588 CurrentA_PointVelocity, CurrentB_PointVelocity,
589 SubStart_PointVelocity, CurrentE_Point,
590 ApproxIntersecPointV, substep_no, substep_no_p, depth);
591 errmsg << G4endl << " * Location: " << MethodName
592 << "- After EndIf(Intersects_AF)" << G4endl;
593 errmsg << " * Bool flags: Recalculated = " << recalculatedB
594 << " Intersects_AF = " << Intersects_AF
595 << " Intersects_FB = " << Intersects_FB << G4endl;
596 errmsg << " * Number of calls to MLL:EIP= " << fNumCalls << G4endl;
597 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
598 }
599 if( restoredFullEndpoint )
600 {
601 fin_section_depth[depth] = restoredFullEndpoint;
602 restoredFullEndpoint = false;
603 }
604 } // EndIf ( E is close enough to the curve, ie point F. )
605 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
606
607#ifdef G4DEBUG_FIELD
608 if( trigger_substepno_print == 0)
609 {
610 trigger_substepno_print= fWarnSteps - 20;
611 }
612
613 if( substep_no >= trigger_substepno_print )
614 {
615 G4cout << "Difficulty in converging in " << MethodName
616 << " Substep no = " << substep_no << G4endl;
617 if( substep_no == trigger_substepno_print )
618 {
619 G4cout << " Start: ";
620 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
621 -1.0, NewSafety, 0 );
622
623 G4cout << " ** Change records: " << G4endl;
624 G4cout << "endPoints A (start) and B (end): combined changes of AB intervals" << G4endl;
625 G4LocatorChangeRecord::ReportEndChanges(G4cout, endChangeA, endChangeB );
626 }
627 G4cout << " Point A: ";
628 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
629 -1.0, NewSafety, substep_no-1 );
630 G4cout << " Point B: ";
631 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
632 -1.0, NewSafety, substep_no );
633 }
634#endif
635 ++substep_no;
636 ++substep_no_p;
637
638 } while ( ( ! found_approximate_intersection )
639 && ( ! there_is_no_intersection )
640 && ( substep_no_p <= param_substeps) ); // UNTIL found or
641 // failed param substep
642
643 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
644 {
645 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
646 - SubStart_PointVelocity.GetCurveLength());
647 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
648 - SubStart_PointVelocity.GetCurveLength());
649
650 G4double distAB = -1;
651 //
652 // Is progress is too slow, and is it possible to go deeper?
653 // If so, then *halve the step*
654 // ==============
655 if( (did_len < fraction_done*all_len)
656 && (depth<max_depth) && (!sub_final_section) )
657 {
658#ifdef G4DEBUG_FIELD
659 static G4ThreadLocal unsigned int numSplits=0; // For debugging only
660 biggest_depth = std::max(depth, biggest_depth);
661 ++numSplits;
662#endif
663 Second_half = false;
664 ++depth;
665 first_section = false;
666
667 G4double Sub_len = (all_len-did_len)/(2.);
668 G4FieldTrack midPoint = CurrentA_PointVelocity;
669 G4bool fullAdvance=
670 integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep);
671
673 if( fullAdvance ) { ++fNumAdvanceFull; }
674
675 G4double lenAchieved=
676 midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength();
677
678 const G4double adequateFraction = (1.0-CLHEP::perThousand);
679 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
680 if ( goodAdvance ) { ++fNumAdvanceGood; }
681
682#ifdef G4DEBUG_FIELD
683 else // !goodAdvance
684 {
685 G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth
686 << " total length achieved = " << lenAchieved << " of "
687 << Sub_len << " fraction= ";
688 if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; }
689 else { G4cout << "DivByZero"; }
690 G4cout << " Good-enough fraction = " << adequateFraction;
691 G4cout << " Number of call to mll = " << fNumCalls
692 << " iteration " << substep_no
693 << " inner = " << substep_no_p << G4endl;
694 G4cout << " Epsilon of integration = " << fiEpsilonStep << G4endl;
695 G4cout << " State at start is = " << CurrentA_PointVelocity
696 << G4endl
697 << " at end (midpoint)= " << midPoint << G4endl;
698 G4cout << " Particle mass = " << midPoint.GetRestMass() << G4endl;
699
700 G4EquationOfMotion *equation = integrDriver->GetEquationOfMotion();
701 ReportFieldValue( CurrentA_PointVelocity, "start", equation );
702 ReportFieldValue( midPoint, "midPoint", equation );
703 G4cout << " Original Start = "
704 << CurveStartPointVelocity << G4endl;
705 G4cout << " Original End = "
706 << CurveEndPointVelocity << G4endl;
707 G4cout << " Original TrialPoint = "
708 << TrialPoint << G4endl;
709 G4cout << " (this is STRICT mode) "
710 << " num Splits= " << numSplits;
711 G4cout << G4endl;
712 }
713#endif
714
715 *ptrInterMedFT[depth] = midPoint;
716 CurrentB_PointVelocity = midPoint;
717
718 if (fCheckMode)
719 {
720 ++eventCount;
721 endChangeB.push_back(
723 substep_no, eventCount, midPoint) );
724 }
725
726 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
727 //
728 SubStart_PointVelocity = CurrentA_PointVelocity;
729
730 // Find new trial intersection point needed at start of the loop
731 //
732 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
733 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
734
735 G4ThreeVector PointGe;
737 G4bool Intersects_AB = IntersectChord(Point_A, Point_B,
738 NewSafety, previousSafety,
739 previousSftOrigin, distAB,
740 PointGe);
741 if( Intersects_AB )
742 {
743 last_AF_intersection = Intersects_AB;
744 CurrentE_Point = PointGe;
745 fin_section_depth[depth] = true;
746
747 validIntersectP = true; // 'E' has been updated.
748
749 G4bool validNormalAB;
750 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
751 validNormalAtE = validNormalAB;
752 }
753 else
754 {
755 // No intersection found for first part of curve
756 // (CurrentA,InterMedPoint[depth]). Go to the second part
757 //
758 Second_half = true;
759
760 validIntersectP= false; // No new 'E' chord intersection found
761 }
762 } // if did_len
763
764 unsigned int levelPops = 0;
765
766 G4bool unfinished = Second_half;
767 while ( unfinished && (depth>0) ) // Loop checking, 07.10.2016, JA
768 {
769 // Second part of curve (InterMed[depth],Intermed[depth-1]))
770 // On the depth-1 level normally we are on the 'second_half'
771
772 ++levelPops;
773
774 // Find new trial intersection point needed at start of the loop
775 //
776 SubStart_PointVelocity = *ptrInterMedFT[depth];
777 CurrentA_PointVelocity = *ptrInterMedFT[depth];
778 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
779
780 if (fCheckMode)
781 {
782 ++eventCount;
784 substep_no, eventCount, CurrentA_PointVelocity);
785 endChangeA.push_back( chngPop_a );
787 substep_no, eventCount, CurrentB_PointVelocity);
788 endChangeB.push_back( chngPop_b );
789 }
790
791 // Ensure that the new endpoints are not further apart in space
792 // than on the curve due to different errors in the integration
793 //
794 G4int errorEndPt = -1;
795 G4bool recalculatedB= false;
796 if( driverReIntegrates )
797 {
798 // Ensure that the new endpoints are not further apart in space
799 // than on the curve due to different errors in the integration
800 //
801 G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity;
802 recalculatedB =
803 CheckAndReEstimateEndpoint( CurrentA_PointVelocity,
804 CurrentB_PointVelocity,
805 RevisedEndPointFT,
806 errorEndPt );
807 if( recalculatedB )
808 {
809 CurrentB_PointVelocity = RevisedEndPointFT; // Use it !
810
811 if ( depth == 1 )
812 {
813 recalculatedEndPoint = true;
814 IntersectedOrRecalculatedFT = RevisedEndPointFT;
815 // So that we can return it, if it is the endpoint!
816 }
817 }
818 else
819 {
820 if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() )
821 errorEndPt = 2;
822 }
823
824 if (fCheckMode)
825 {
826 ++eventCount;
827 endChangeB.push_back(
829 substep_no, eventCount, RevisedEndPointFT));
830 }
831 }
832 if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
833 {
834 std::ostringstream errmsg;
835 ReportReversedPoints(errmsg,
836 CurveStartPointVelocity, CurveEndPointVelocity,
837 NewSafety, fiEpsilonStep,
838 CurrentA_PointVelocity, CurrentB_PointVelocity,
839 SubStart_PointVelocity, CurrentE_Point,
840 ApproxIntersecPointV, substep_no, substep_no_p, depth);
841 errmsg << " * Location: " << MethodName << "- Second-Half" << G4endl;
842 errmsg << " * Recalculated = " << recalculatedB << G4endl; // false
843 G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
844 }
845
846 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
847 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
848 G4ThreeVector PointGi;
850 G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety,
851 previousSafety,
852 previousSftOrigin, distAB,
853 PointGi);
854 if( Intersects_AB )
855 {
856 last_AF_intersection = Intersects_AB;
857 CurrentE_Point = PointGi;
858
859 validIntersectP = true; // 'E' has been updated.
860 NormalAtEntry = GetSurfaceNormal( PointGi, validNormalAtE );
861 }
862 else
863 {
864 validIntersectP = false; // No new 'E' chord intersection found
865 if( depth == 1)
866 {
867 there_is_no_intersection = true;
868 }
869 }
870 depth--;
871 fin_section_depth[depth] = true;
872 unfinished = !validIntersectP;
873 }
874#ifdef G4DEBUG_FIELD
875 if( ! ( validIntersectP || there_is_no_intersection ) )
876 {
877 // What happened ??
878 G4cout << "MLL - WARNING Potential FAILURE: Conditions not met!"
879 << G4endl
880 << " Depth = " << depth << G4endl
881 << " Levels popped = " << levelPops
882 << " Num Substeps= " << substep_no << G4endl;
883 G4cout << " Found intersection= " << found_approximate_intersection
884 << G4endl;
885 G4cout << " Progress report: -- " << G4endl;
887 CurveStartPointVelocity, CurveEndPointVelocity,
888 substep_no, CurrentA_PointVelocity,
889 CurrentB_PointVelocity,
890 NewSafety, depth );
891 G4cout << G4endl;
892 }
893#endif
894 } // if(!found_aproximate_intersection)
895
896 assert( validIntersectP || there_is_no_intersection
897 || found_approximate_intersection);
898
899 } while ( ( ! found_approximate_intersection )
900 && ( ! there_is_no_intersection )
901 && ( substep_no <= fMaxSteps) ); // UNTIL found or failed
902
903 if( substep_no > max_no_seen )
904 {
905 max_no_seen = substep_no;
906#ifdef G4DEBUG_FIELD
907 if( max_no_seen > fWarnSteps )
908 {
909 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
910 }
911#endif
912 }
913
914 if( !there_is_no_intersection && !found_approximate_intersection )
915 {
916 if( substep_no >= fMaxSteps)
917 {
918 // Since we cannot go further (yet), we return as far as we have gone
919
920 recalculatedEndPoint = true;
921 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
922 found_approximate_intersection = false;
923
924 std::ostringstream message;
925 message << G4endl;
926 message << "Convergence is requiring too many substeps: "
927 << substep_no << " ( limit = "<< fMaxSteps << ")"
928 << G4endl
929 << " Abandoning effort to intersect. " << G4endl << G4endl;
930 message << " Number of calls to MLL: " << fNumCalls;
931 message << " iteration = " << substep_no <<G4endl << G4endl;
932
933 message.precision( 12 );
934 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
935 G4double full_len = CurveEndPointVelocity.GetCurveLength();
936 message << " Undertaken only length: " << done_len
937 << " out of " << full_len << " required." << G4endl
938 << " Remaining length = " << full_len - done_len;
939
940 message << " Start and end-point of requested Step:" << G4endl;
941 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
942 -1.0, NewSafety, 0, message, -1 );
943 message << " Start and end-point of current Sub-Step:" << G4endl;
944 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
945 -1.0, NewSafety, substep_no-1, message, -1 );
946 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
947 -1.0, NewSafety, substep_no, message, -1 );
948
949 G4Exception(MethodName, "GeomNav0003", JustWarning, message);
950 }
951 else if( substep_no >= fWarnSteps)
952 {
953 std::ostringstream message;
954 message << "Many substeps while trying to locate intersection."
955 << G4endl
956 << " Undertaken length: "
957 << CurrentB_PointVelocity.GetCurveLength()
958 << " - Needed: " << substep_no << " substeps." << G4endl
959 << " Warning number = " << fWarnSteps
960 << " and maximum substeps = " << fMaxSteps;
961 G4Exception(MethodName, "GeomNav1002", JustWarning, message);
962 }
963 }
964
965 return (!there_is_no_intersection) && found_approximate_intersection;
966 // Success or failure
967}
@ FatalException
CLHEP::Hep3Vector G4ThreeVector
int G4int
Definition: G4Types.hh:85
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
double mag2() const
double mag() const
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
G4VIntegrationDriver * GetIntegrationDriver()
G4double GetRestMass() const
static std::ostream & ReportEndChanges(std::ostream &os, const G4LocatorChangeLogger &startA, const G4LocatorChangeLogger &endB)
static std::ostream & ReportEndChanges(std::ostream &os, const std::vector< G4LocatorChangeRecord > &startA, const std::vector< G4LocatorChangeRecord > &endB)
void ReportFieldValue(const G4FieldTrack &locationPV, const char *nameLoc, const G4EquationOfMotion *equation)
unsigned long int fNumCalls
unsigned long int fNumAdvanceTrials
unsigned long int fNumAdvanceGood
unsigned long int fNumAdvanceFull
virtual G4bool DoesReIntegrate() const =0
G4ChordFinder * GetChordFinderFor()
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
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)
G4double GetEpsilonStepFor()
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, G4double tolerance, unsigned long int numCalls)
G4bool GetAdjustementOfFoundIntersection()
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
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)
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)
static constexpr double mm
Definition: SystemOfUnits.h:96
static constexpr double perThousand
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T sqr(const T &x)
Definition: templates.hh:128
#define G4ThreadLocal
Definition: tls.hh:77

References G4VIntersectionLocator::AdjustmentOfFoundIntersection(), G4ChordFinder::ApproxCurvePointV(), G4VIntersectionLocator::CheckAndReEstimateEndpoint(), G4VIntegrationDriver::DoesReIntegrate(), CLHEP::Hep3Vector::dot(), FatalException, G4VIntersectionLocator::fCheckMode, G4VIntersectionLocator::fiDeltaIntersection, G4VIntersectionLocator::fiEpsilonStep, fMaxSteps, fNumAdvanceFull, fNumAdvanceGood, fNumAdvanceTrials, fNumCalls, G4VIntersectionLocator::fVerboseLevel, fWarnSteps, G4cerr, G4cout, G4endl, G4Exception(), G4ThreadLocal, G4VIntersectionLocator::GetAdjustementOfFoundIntersection(), G4VIntersectionLocator::GetChordFinderFor(), G4FieldTrack::GetCurveLength(), G4VIntersectionLocator::GetEpsilonStepFor(), G4ChordFinder::GetIntegrationDriver(), G4FieldTrack::GetMomentumDir(), G4FieldTrack::GetMomentumDirection(), G4VIntersectionLocator::GetNavigatorFor(), G4FieldTrack::GetPosition(), G4FieldTrack::GetRestMass(), G4VIntersectionLocator::GetSurfaceNormal(), G4VIntersectionLocator::IntersectChord(), JustWarning, G4VIntersectionLocator::kCarTolerance, G4LocatorChangeRecord::kInitialisingCL, G4LocatorChangeRecord::kInsertingMidPoint, G4LocatorChangeRecord::kIntersectsAF, G4LocatorChangeRecord::kIntersectsFB, G4LocatorChangeRecord::kInvalidCL, G4LocatorChangeRecord::kLevelPop, G4LocatorChangeRecord::kNoIntersectAForFB, G4LocatorChangeRecord::kRecalculatedB, G4LocatorChangeRecord::kRecalculatedBagn, G4Navigator::LocateGlobalPointWithinVolume(), CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), G4INCL::Math::max(), max_depth, CLHEP::mm, CLHEP::perThousand, G4VIntersectionLocator::printStatus(), ptrInterMedFT, G4LocatorChangeLogger::ReportEndChanges(), G4LocatorChangeRecord::ReportEndChanges(), ReportFieldValue(), G4VIntersectionLocator::ReportImmediateHit(), G4VIntersectionLocator::ReportProgress(), G4VIntersectionLocator::ReportReversedPoints(), G4VIntersectionLocator::ReportTrialStep(), G4FieldTrack::SetMomentum(), G4FieldTrack::SetPosition(), and sqr().

◆ GetAdjustementOfFoundIntersection()

G4bool G4VIntersectionLocator::GetAdjustementOfFoundIntersection ( )
inlineinherited

◆ GetCheckMode()

G4bool G4VIntersectionLocator::GetCheckMode ( )
inlineinherited

◆ GetChordFinderFor()

G4ChordFinder * G4VIntersectionLocator::GetChordFinderFor ( )
inlineinherited

◆ GetDeltaIntersectionFor()

G4double G4VIntersectionLocator::GetDeltaIntersectionFor ( )
inlineinherited

◆ GetEpsilonStepFor()

G4double G4VIntersectionLocator::GetEpsilonStepFor ( )
inlineinherited

◆ GetGlobalSurfaceNormal()

G4ThreeVector G4VIntersectionLocator::GetGlobalSurfaceNormal ( const G4ThreeVector CurrentE_Point,
G4bool validNormal 
)
protectedinherited

Definition at line 563 of file G4VIntersectionLocator.cc.

566{
567 G4ThreeVector localNormal = GetLocalSurfaceNormal(CurrentE_Point,validNormal);
568 G4AffineTransform localToGlobal = // Must use the same Navigator !!
570 G4ThreeVector globalNormal = localToGlobal.TransformAxis( localNormal );
571
572#ifdef G4DEBUG_FIELD
573 if( validNormal && ( std::fabs(globalNormal.mag2() - 1.0) > perThousand ) )
574 {
575 std::ostringstream message;
576 message << "**************************************************************"
577 << G4endl;
578 message << " Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
579 << G4endl;
580 message << " * Constituents: " << G4endl;
581 message << " Local Normal= " << localNormal << G4endl;
582 message << " Transform: " << G4endl
583 << " Net Translation= " << localToGlobal.NetTranslation()
584 << G4endl
585 << " Net Rotation = " << localToGlobal.NetRotation()
586 << G4endl;
587 message << " * Result: " << G4endl;
588 message << " Global Normal= " << localNormal << G4endl;
589 message << "**************************************************************";
590 G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
591 "GeomNav1002", JustWarning, message);
592 }
593#endif
594
595 return globalNormal;
596}
static constexpr double perThousand
Definition: G4SIunits.hh:326
G4ThreeVector NetTranslation() const
G4RotationMatrix NetRotation() const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
const G4AffineTransform GetLocalToGlobalTransform() const
G4ThreeVector GetLocalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)

References G4VIntersectionLocator::fHelpingNavigator, G4endl, G4Exception(), G4VIntersectionLocator::GetLocalSurfaceNormal(), G4Navigator::GetLocalToGlobalTransform(), JustWarning, CLHEP::Hep3Vector::mag2(), G4AffineTransform::NetRotation(), G4AffineTransform::NetTranslation(), perThousand, and G4AffineTransform::TransformAxis().

Referenced by G4VIntersectionLocator::AdjustmentOfFoundIntersection().

◆ GetLastSurfaceNormal()

G4ThreeVector G4VIntersectionLocator::GetLastSurfaceNormal ( const G4ThreeVector intersectPoint,
G4bool validNormal 
) const
privateinherited

Definition at line 602 of file G4VIntersectionLocator.cc.

605{
606 G4ThreeVector normalVec;
607 G4bool validNorm;
608 normalVec = fiNavigator->GetGlobalExitNormal( intersectPoint, &validNorm );
609 normalIsValid = validNorm;
610
611 return normalVec;
612}
virtual G4ThreeVector GetGlobalExitNormal(const G4ThreeVector &point, G4bool *valid)

References G4VIntersectionLocator::fiNavigator, and G4Navigator::GetGlobalExitNormal().

Referenced by G4VIntersectionLocator::GetSurfaceNormal().

◆ GetLocalSurfaceNormal()

G4ThreeVector G4VIntersectionLocator::GetLocalSurfaceNormal ( const G4ThreeVector CurrentE_Point,
G4bool validNormal 
)
privateinherited

Definition at line 384 of file G4VIntersectionLocator.cc.

386{
387 G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0.0));
388 G4VPhysicalVolume* located;
389
390 validNormal = false;
392 located = fHelpingNavigator->LocateGlobalPointAndSetup( CurrentE_Point );
393
394 delete fpTouchable;
396
397 // To check if we can use GetGlobalExitNormal()
398 //
399 G4ThreeVector localPosition = fpTouchable->GetHistory()
400 ->GetTopTransform().TransformPoint(CurrentE_Point);
401
402 // Issue: in the case of coincident surfaces, this version does not recognise
403 // which side you are located onto (can return vector with wrong sign.)
404 // TO-DO: use direction (of chord) to identify volume we will be "entering"
405
406 if( located != 0)
407 {
408 G4LogicalVolume* pLogical= located->GetLogicalVolume();
409 G4VSolid* pSolid;
410
411 if( (pLogical != nullptr) && ( (pSolid=pLogical->GetSolid()) != nullptr ) )
412 {
413 if ( ( pSolid->Inside(localPosition)==kSurface )
414 || ( pSolid->DistanceToOut(localPosition) < 1000.0 * kCarTolerance ) )
415 {
416 Normal = pSolid->SurfaceNormal(localPosition);
417 validNormal = true;
418
419#ifdef G4DEBUG_FIELD
420 if( std::fabs(Normal.mag2() - 1.0 ) > CLHEP::perThousand)
421 {
422 G4cerr << "PROBLEM in G4VIntersectionLocator::GetLocalSurfaceNormal."
423 << G4endl;
424 G4cerr << " Normal is not unit - mag=" << Normal.mag() << G4endl;
425 G4cerr << " at trial local point " << CurrentE_Point << G4endl;
426 G4cerr << " Solid is " << *pSolid << G4endl;
427 }
428#endif
429 }
430 }
431 }
432 return Normal;
433}
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4VSolid * GetSolid() const
const G4AffineTransform & GetTopTransform() const
G4TouchableHistory * CreateTouchableHistory() const
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=nullptr, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:132
void SetWorldVolume(G4VPhysicalVolume *pWorld)
const G4NavigationHistory * GetHistory() const
G4TouchableHistory * fpTouchable
G4LogicalVolume * GetLogicalVolume() const
virtual EInside Inside(const G4ThreeVector &p) const =0
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const =0
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
@ kSurface
Definition: geomdefs.hh:69

References G4Navigator::CreateTouchableHistory(), G4VSolid::DistanceToOut(), G4VIntersectionLocator::fHelpingNavigator, G4VIntersectionLocator::fpTouchable, G4cerr, G4endl, G4TouchableHistory::GetHistory(), G4VPhysicalVolume::GetLogicalVolume(), G4VIntersectionLocator::GetNavigatorFor(), G4LogicalVolume::GetSolid(), G4NavigationHistory::GetTopTransform(), G4VSolid::Inside(), G4VIntersectionLocator::kCarTolerance, kSurface, G4Navigator::LocateGlobalPointAndSetup(), CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), CLHEP::perThousand, G4Navigator::SetWorldVolume(), G4VSolid::SurfaceNormal(), and G4AffineTransform::TransformPoint().

Referenced by G4VIntersectionLocator::GetGlobalSurfaceNormal().

◆ GetNavigatorFor()

G4Navigator * G4VIntersectionLocator::GetNavigatorFor ( )
inlineinherited

◆ GetSurfaceNormal()

G4ThreeVector G4VIntersectionLocator::GetSurfaceNormal ( const G4ThreeVector CurrentInt_Point,
G4bool validNormal 
)
protectedinherited

Definition at line 518 of file G4VIntersectionLocator.cc.

521{
522 G4ThreeVector NormalAtEntry; // ( -10. , -10., -10. );
523
524 G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
525 G4bool validNormalLast;
526
527 // Relies on a call to Navigator::ComputeStep in IntersectChord before
528 // this call
529 //
530 NormalAtEntryLast = GetLastSurfaceNormal( CurrentInt_Point, validNormalLast );
531 // May return valid=false in cases, including
532 // - if the candidate volume was not found (eg exiting world), or
533 // - a replica was involved -- determined the step size.
534 // (This list is not complete.)
535
536#ifdef G4DEBUG_FIELD
537 if ( validNormalLast
538 && ( std::fabs(NormalAtEntryLast.mag2() - 1.0) > perThousand ) )
539 {
540 std::ostringstream message;
541 message << "PROBLEM: Normal is not unit - magnitude = "
542 << NormalAtEntryLast.mag() << G4endl;
543 message << " at trial intersection point " << CurrentInt_Point << G4endl;
544 message << " Obtained from Get *Last* Surface Normal.";
545 G4Exception("G4VIntersectionLocator::GetSurfaceNormal()",
546 "GeomNav1002", JustWarning, message);
547 }
548#endif
549
550 if( validNormalLast )
551 {
552 NormalAtEntry = NormalAtEntryLast;
553 }
554 validNormal = validNormalLast;
555
556 return NormalAtEntry;
557}
G4ThreeVector GetLastSurfaceNormal(const G4ThreeVector &intersectPoint, G4bool &validNormal) const

References G4endl, G4Exception(), G4VIntersectionLocator::GetLastSurfaceNormal(), JustWarning, CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), and perThousand.

Referenced by G4BrentLocator::EstimateIntersectionPoint(), EstimateIntersectionPoint(), and G4SimpleLocator::EstimateIntersectionPoint().

◆ GetVerboseFor()

G4int G4VIntersectionLocator::GetVerboseFor ( )
inlineinherited

◆ IntersectChord()

G4bool G4VIntersectionLocator::IntersectChord ( const G4ThreeVector StartPointA,
const G4ThreeVector EndPointB,
G4double NewSafety,
G4double PreviousSafety,
G4ThreeVector PreviousSftOrigin,
G4double LinearStepLength,
G4ThreeVector IntersectionPoint,
G4bool calledNavigator = nullptr 
)
inlineinherited

◆ LocateGlobalPointWithinVolumeAndCheck()

G4bool G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck ( const G4ThreeVector pos)
protectedinherited

Definition at line 674 of file G4VIntersectionLocator.cc.

676{
677 G4bool good = true;
679 const G4String
680 MethodName("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()");
681
682 if( fCheckMode )
683 {
684 G4bool navCheck= nav->IsCheckModeActive(); // Recover original value
685 nav->CheckMode(true);
686
687 // Identify the current volume
688
690 G4VPhysicalVolume* motherPhys = startTH->GetVolume();
691 G4VSolid* motherSolid = startTH->GetSolid();
693 G4int motherCopyNo = motherPhys->GetCopyNo();
694
695 // Let's check that the point is inside the current solid
696 G4ThreeVector localPosition = transform.TransformPoint(position);
697 EInside inMother = motherSolid->Inside( localPosition );
698 if( inMother != kInside )
699 {
700 std::ostringstream message;
701 message << "Position located "
702 << ( inMother == kSurface ? " on Surface " : " outside " )
703 << "expected volume" << G4endl
704 << " Safety (from Outside) = "
705 << motherSolid->DistanceToIn(localPosition);
706 G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
707 "GeomNav1002", JustWarning, message);
708 }
709
710 // 1. Simple next step - quick relocation and check result.
711 // nav->LocateGlobalPointWithinVolume( position );
712
713 // 2. Full relocation - to cross-check answer !
715 if( (nextPhysical != motherPhys)
716 || (nextPhysical->GetCopyNo() != motherCopyNo )
717 )
718 {
719 G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
720 "GeomNav1002", JustWarning,
721 "Position located outside expected volume.");
722 }
723 nav->CheckMode(navCheck); // Recover original value
724 }
725 else
726 {
728 }
729 return good;
730}
void CheckMode(G4bool mode)
virtual G4TouchableHistoryHandle CreateTouchableHistoryHandle() const
G4bool IsCheckModeActive() const
const G4AffineTransform & GetGlobalToLocalTransform() const
virtual G4int GetCopyNo() const =0
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
EInside
Definition: geomdefs.hh:67
@ kInside
Definition: geomdefs.hh:70
G4bool transform(G4String &input, const G4String &type)

References G4Navigator::CheckMode(), G4Navigator::CreateTouchableHistoryHandle(), G4VSolid::DistanceToIn(), G4VIntersectionLocator::fCheckMode, G4endl, G4Exception(), G4VPhysicalVolume::GetCopyNo(), G4Navigator::GetGlobalToLocalTransform(), G4VIntersectionLocator::GetNavigatorFor(), G4VSolid::Inside(), G4Navigator::IsCheckModeActive(), JustWarning, kInside, kSurface, G4Navigator::LocateGlobalPointAndSetup(), G4Navigator::LocateGlobalPointWithinVolume(), and G4coutFormatters::anonymous_namespace{G4coutFormatters.cc}::transform().

Referenced by G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport().

◆ LocateGlobalPointWithinVolumeCheckAndReport()

void G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport ( const G4ThreeVector pos,
const G4String CodeLocationInfo,
G4int  CheckMode 
)
protectedinherited

Definition at line 736 of file G4VIntersectionLocator.cc.

740{
741 // Save value of Check mode first
742 G4bool oldCheck = GetCheckMode();
743
745 if( !ok )
746 {
747 std::ostringstream message;
748 message << "Failed point location." << G4endl
749 << " Code Location info: " << CodeLocationInfo;
750 G4Exception("G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport()",
751 "GeomNav1002", JustWarning, message);
752 }
753
754 SetCheckMode( oldCheck );
755}
G4bool LocateGlobalPointWithinVolumeAndCheck(const G4ThreeVector &pos)

References G4endl, G4Exception(), G4VIntersectionLocator::GetCheckMode(), JustWarning, G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck(), and G4VIntersectionLocator::SetCheckMode().

◆ printStatus() [1/2]

void G4VIntersectionLocator::printStatus ( const G4FieldTrack startFT,
const G4FieldTrack currentFT,
G4double  requestStep,
G4double  safety,
G4int  stepNum 
)
inherited

◆ printStatus() [2/2]

void G4VIntersectionLocator::printStatus ( const G4FieldTrack startFT,
const G4FieldTrack currentFT,
G4double  requestStep,
G4double  safety,
G4int  stepNum,
std::ostream &  oss,
G4int  verboseLevel 
)
staticinherited

Definition at line 91 of file G4VIntersectionLocator.cc.

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}
double z() const
double x() const
double y() const

References G4endl, G4FieldTrack::GetCurveLength(), G4FieldTrack::GetMomentum(), G4FieldTrack::GetMomentumDir(), G4FieldTrack::GetPosition(), CLHEP::Hep3Vector::mag(), G4VIntersectionLocator::printStatus(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

◆ ReEstimateEndpoint()

G4FieldTrack G4VIntersectionLocator::ReEstimateEndpoint ( const G4FieldTrack CurrentStateA,
const G4FieldTrack EstimtdEndStateB,
G4double  linearDistSq,
G4double  curveDist 
)
protectedinherited

Definition at line 188 of file G4VIntersectionLocator.cc.

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,
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/*
300#else
301 // Statistics on the RMS value of the corrections
302
303 static G4ThreadLocal G4int noCorrections = 0;
304 ++noCorrections;
305 if( goodAdvance )
306 {
307 static G4ThreadLocal G4double sumCorrectionsSq;
308 sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
309 newEndPoint.GetPosition()).mag2();
310 }
311*/
312#endif
313
314 return retEndPoint;
315}

References G4VIntersectionLocator::fVerboseLevel, G4cout, G4endl, G4Exception(), G4VIntersectionLocator::GetChordFinderFor(), G4FieldTrack::GetCurveLength(), G4VIntersectionLocator::GetEpsilonStepFor(), G4ChordFinder::GetIntegrationDriver(), G4FieldTrack::GetPosition(), JustWarning, G4VIntersectionLocator::kCarTolerance, and CLHEP::Hep3Vector::mag().

Referenced by G4VIntersectionLocator::CheckAndReEstimateEndpoint(), G4BrentLocator::EstimateIntersectionPoint(), and G4SimpleLocator::EstimateIntersectionPoint().

◆ ReportFieldValue()

void G4MultiLevelLocator::ReportFieldValue ( const G4FieldTrack locationPV,
const char *  nameLoc,
const G4EquationOfMotion equation 
)
private

Definition at line 980 of file G4MultiLevelLocator.cc.

983{
984 enum { maxNumFieldComp = 24 };
985
986 G4ThreeVector position = locationPV.GetPosition();
987 G4double startPoint[4] = { position.x(), position.y(), position.z(),
988 locationPV.GetLabTimeOfFlight() };
989 G4double FieldVec[maxNumFieldComp]; // 24 ;
990 for (auto i=0; i<maxNumFieldComp; ++i )
991 {
992 FieldVec[i] = 0.0;
993 }
994 equation->GetFieldValue( startPoint, FieldVec);
995 G4cout << " B-field value (" << nameLoc << ")= "
996 << FieldVec[0] << " " << FieldVec[1] << " " << FieldVec[2];
997 G4double Emag2= G4ThreeVector( FieldVec[3],
998 FieldVec[4],
999 FieldVec[5] ).mag2();
1000 if( Emag2 > 0.0 )
1001 {
1002 G4cout << " Electric = " << FieldVec[3] << " "
1003 << FieldVec[4] << " "
1004 << FieldVec[5]<< G4endl;
1005 }
1006 return;
1007}
void GetFieldValue(const G4double Point[4], G4double Field[]) const
G4double GetLabTimeOfFlight() const

References G4cout, G4endl, G4EquationOfMotion::GetFieldValue(), G4FieldTrack::GetLabTimeOfFlight(), G4FieldTrack::GetPosition(), and CLHEP::Hep3Vector::mag2().

Referenced by EstimateIntersectionPoint().

◆ ReportImmediateHit()

void G4VIntersectionLocator::ReportImmediateHit ( const char *  MethodName,
const G4ThreeVector StartPosition,
const G4ThreeVector TrialPoint,
G4double  tolerance,
unsigned long int  numCalls 
)
protectedinherited

Definition at line 845 of file G4VIntersectionLocator.cc.

850{
851 static G4ThreadLocal unsigned int occurredOnTop= 0;
852 static G4ThreadLocal G4ThreeVector* ptrLast = nullptr;
853 if( ptrLast == nullptr )
854 {
855 ptrLast= new G4ThreeVector( DBL_MAX, DBL_MAX, DBL_MAX );
856 G4AutoDelete::Register(ptrLast);
857 }
858 G4ThreeVector &lastStart= *ptrLast;
859
860 if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
861 {
862 static G4ThreadLocal unsigned int numUnmoved = 0;
863 static G4ThreadLocal unsigned int numStill = 0; // Still at same point
864
865 G4cout << "Intersection F == start A in " << MethodName;
866 G4cout << "Start Point: " << StartPosition << G4endl;
867 G4cout << " Start-Trial: " << TrialPoint - StartPosition;
868 G4cout << " Start-last: " << StartPosition - lastStart;
869
870 if( (StartPosition - lastStart).mag() < tolerance )
871 {
872 // We are at position of last 'Start' position - ie unmoved
873 ++numUnmoved;
874 ++numStill;
875 G4cout << " { Unmoved: " << " still#= " << numStill
876 << " total # = " << numUnmoved << " } - ";
877 }
878 else
879 {
880 numStill = 0;
881 }
882 G4cout << " Occurred: " << ++occurredOnTop;
883 G4cout << " out of total calls= " << numCalls;
884 G4cout << G4endl;
885 lastStart = StartPosition;
886 }
887} // End of ReportImmediateHit()
void Register(T *inst)
Definition: G4AutoDelete.hh:65
#define DBL_MAX
Definition: templates.hh:62

References DBL_MAX, G4cout, G4endl, G4ThreadLocal, and G4AutoDelete::Register().

Referenced by EstimateIntersectionPoint().

◆ ReportProgress()

void G4VIntersectionLocator::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 
)
protectedinherited

Definition at line 813 of file G4VIntersectionLocator.cc.

822{
823 oss << "ReportProgress: Current status of intersection search: " << G4endl;
824 if( depth > 0 ) oss << " Depth= " << depth;
825 oss << " Substep no = " << substep_no << G4endl;
826 G4int verboseLevel = 5;
827 G4double safetyPrev = -1.0; // Add as argument ?
828
829 printStatus( StartPointVel, EndPointVel, -1.0, -1.0, -1,
830 oss, verboseLevel);
831 oss << " * Start and end-point of requested Step:" << G4endl;
832 oss << " ** State of point A: ";
833 printStatus( A_PtVel, A_PtVel, -1.0, safetyPrev, substep_no-1,
834 oss, verboseLevel);
835 oss << " ** State of point B: ";
836 printStatus( A_PtVel, B_PtVel, -1.0, safetyLast, substep_no,
837 oss, verboseLevel);
838}

References G4endl, and G4VIntersectionLocator::printStatus().

Referenced by EstimateIntersectionPoint().

◆ ReportReversedPoints()

void G4VIntersectionLocator::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 
)
protectedinherited

Definition at line 761 of file G4VIntersectionLocator.cc.

772{
773 // Expect that 'msg' can hold the name of the calling method
774
775 // FieldTrack 'points' A and B have been tangled
776 // Whereas A should be before B, it is found that curveLen(B) < curveLen(A)
777 G4int verboseLevel= 5;
778 G4double curveDist = B_PtVel.GetCurveLength() - A_PtVel.GetCurveLength();
779 G4VIntersectionLocator::printStatus( A_PtVel, B_PtVel,
780 -1.0, NewSafety, substep_no, msg, verboseLevel );
781 msg << "Error in advancing propagation." << G4endl
782 << " The final curve point is NOT further along"
783 << " than the original!" << G4endl
784 << " Going *backwards* from len(A) = " << A_PtVel.GetCurveLength()
785 << " to len(B) = " << B_PtVel.GetCurveLength() << G4endl
786 << " Curve distance is " << curveDist / CLHEP::millimeter << " mm "
787 << G4endl
788 << " Point A' (start) is " << A_PtVel << G4endl
789 << " Point B' (end) is " << B_PtVel << G4endl;
790 msg << " fEpsStep= " << epsStep << G4endl << G4endl;
791
792 G4int oldprc = msg.precision(20);
793 msg << " In full precision, the position, momentum, E_kin, length, rest mass "
794 << " ... are: " << G4endl;
795 msg << " Point A[0] (Curve start) is " << StartPointVel << G4endl
796 << " Point S (Sub start) is " << SubStart_PtVel
797 << " Point A' (Current start) is " << A_PtVel << G4endl
798 << " Point E (Trial Point) is " << E_Point << G4endl
799 << " Point F (Intersection) is " << ApproxIntersecPointV << G4endl
800 << " Point B' (Current end) is " << B_PtVel << G4endl
801 << " Point B[0] (Curve end) is " << EndPointVel << G4endl
802 << G4endl
803 << " LocateIntersection parameters are : " << G4endl
804 << " Substep no (total) = " << substep_no << G4endl
805 << " Substep no = " << substep_no_p << " at depth= " << depth;
806 msg.precision(oldprc);
807}
static constexpr double millimeter
Definition: SystemOfUnits.h:63

References G4endl, G4FieldTrack::GetCurveLength(), CLHEP::millimeter, and G4VIntersectionLocator::printStatus().

Referenced by EstimateIntersectionPoint().

◆ ReportStatistics()

void G4MultiLevelLocator::ReportStatistics ( )

Definition at line 969 of file G4MultiLevelLocator.cc.

970{
971 G4cout << " Number of calls = " << fNumCalls << G4endl;
972 G4cout << " Number of split level ('advances'): "
974 G4cout << " Number of full advances: "
976 G4cout << " Number of good advances: "
978}

References fNumAdvanceFull, fNumAdvanceGood, fNumAdvanceTrials, fNumCalls, G4cout, and G4endl.

Referenced by ~G4MultiLevelLocator().

◆ ReportTrialStep()

void G4VIntersectionLocator::ReportTrialStep ( G4int  step_no,
const G4ThreeVector ChordAB_v,
const G4ThreeVector ChordEF_v,
const G4ThreeVector NewMomentumDir,
const G4ThreeVector NormalAtEntry,
G4bool  validNormal 
)
protectedinherited

Definition at line 618 of file G4VIntersectionLocator.cc.

624{
625 G4double ABchord_length = ChordAB_v.mag();
626 G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry );
627 G4double MomDir_dot_ABchord;
628 MomDir_dot_ABchord = (1.0 / ABchord_length) * NewMomentumDir.dot( ChordAB_v );
629
630 std::ostringstream outStream;
631 outStream << std::setw(6) << " Step# "
632 << std::setw(17) << " |ChordEF|(mag)" << " "
633 << std::setw(18) << " uMomentum.Normal" << " "
634 << std::setw(18) << " uMomentum.ABdir " << " "
635 << std::setw(16) << " AB-dist " << " "
636 << " Chord Vector (EF) "
637 << G4endl;
638 outStream.precision(7);
639 outStream << " " << std::setw(5) << step_no
640 << " " << std::setw(18) << ChordEF_v.mag()
641 << " " << std::setw(18) << MomDir_dot_Norm
642 << " " << std::setw(18) << MomDir_dot_ABchord
643 << " " << std::setw(12) << ABchord_length
644 << " " << ChordEF_v
645 << G4endl;
646 outStream << " MomentumDir= " << " " << NewMomentumDir
647 << " Normal at Entry E= " << NormalAtEntry
648 << " AB chord = " << ChordAB_v
649 << G4endl;
650 G4cout << outStream.str();
651
652 if( ( std::fabs(NormalAtEntry.mag2() - 1.0) > perThousand ) )
653 {
654 std::ostringstream message;
655 message << "Normal is not unit - mag= " << NormalAtEntry.mag() << G4endl
656 << " ValidNormalAtE = " << validNormal;
657 G4Exception("G4VIntersectionLocator::ReportTrialStep()",
658 "GeomNav1002", JustWarning, message);
659 }
660 return;
661}

References CLHEP::Hep3Vector::dot(), G4cout, G4endl, G4Exception(), JustWarning, CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), and perThousand.

Referenced by G4BrentLocator::EstimateIntersectionPoint(), EstimateIntersectionPoint(), and G4SimpleLocator::EstimateIntersectionPoint().

◆ SetCheckMode()

void G4VIntersectionLocator::SetCheckMode ( G4bool  value)
inlineinherited

◆ SetChordFinderFor()

void G4VIntersectionLocator::SetChordFinderFor ( G4ChordFinder fCFinder)
inlineinherited

◆ SetDeltaIntersectionFor()

void G4VIntersectionLocator::SetDeltaIntersectionFor ( G4double  deltaIntersection)
inlineinherited

◆ SetEpsilonStepFor()

void G4VIntersectionLocator::SetEpsilonStepFor ( G4double  EpsilonStep)
inlineinherited

◆ SetMaxSteps()

void G4MultiLevelLocator::SetMaxSteps ( unsigned int  valMax)
inline

Definition at line 71 of file G4MultiLevelLocator.hh.

71{ fMaxSteps= valMax; }

References fMaxSteps.

Referenced by G4MultiLevelLocator().

◆ SetNavigatorFor()

void G4VIntersectionLocator::SetNavigatorFor ( G4Navigator fNavigator)
inlineinherited

◆ SetSafetyParametersFor()

void G4VIntersectionLocator::SetSafetyParametersFor ( G4bool  UseSafety)
inlineinherited

◆ SetVerboseFor()

void G4VIntersectionLocator::SetVerboseFor ( G4int  fVerbose)
inlineinherited

◆ SetWarnSteps()

void G4MultiLevelLocator::SetWarnSteps ( unsigned int  valWarn)
inline

Definition at line 72 of file G4MultiLevelLocator.hh.

72{ fWarnSteps= valWarn; }

References fWarnSteps.

Referenced by G4MultiLevelLocator().

Field Documentation

◆ fCheckMode

G4bool G4VIntersectionLocator::fCheckMode = false
protectedinherited

◆ fHelpingNavigator

G4Navigator* G4VIntersectionLocator::fHelpingNavigator
protectedinherited

◆ fiChordFinder

G4ChordFinder* G4VIntersectionLocator::fiChordFinder = nullptr
protectedinherited

Definition at line 261 of file G4VIntersectionLocator.hh.

◆ fiDeltaIntersection

G4double G4VIntersectionLocator::fiDeltaIntersection = -1.0
protectedinherited

◆ fiEpsilonStep

G4double G4VIntersectionLocator::fiEpsilonStep = -1.0
protectedinherited

◆ fiNavigator

G4Navigator* G4VIntersectionLocator::fiNavigator
protectedinherited

◆ fiUseSafety

G4bool G4VIntersectionLocator::fiUseSafety = false
protectedinherited

Definition at line 257 of file G4VIntersectionLocator.hh.

◆ fMaxSteps

unsigned int G4MultiLevelLocator::fMaxSteps = 10000
private

Definition at line 83 of file G4MultiLevelLocator.hh.

Referenced by EstimateIntersectionPoint(), and SetMaxSteps().

◆ fNumAdvanceFull

unsigned long int G4MultiLevelLocator::fNumAdvanceFull = 0
private

Definition at line 92 of file G4MultiLevelLocator.hh.

Referenced by EstimateIntersectionPoint(), and ReportStatistics().

◆ fNumAdvanceGood

unsigned long int G4MultiLevelLocator::fNumAdvanceGood = 0
private

Definition at line 93 of file G4MultiLevelLocator.hh.

Referenced by EstimateIntersectionPoint(), and ReportStatistics().

◆ fNumAdvanceTrials

unsigned long int G4MultiLevelLocator::fNumAdvanceTrials = 0
private

Definition at line 94 of file G4MultiLevelLocator.hh.

Referenced by EstimateIntersectionPoint(), and ReportStatistics().

◆ fNumCalls

unsigned long int G4MultiLevelLocator::fNumCalls = 0
private

Definition at line 91 of file G4MultiLevelLocator.hh.

Referenced by EstimateIntersectionPoint(), and ReportStatistics().

◆ fpTouchable

G4TouchableHistory* G4VIntersectionLocator::fpTouchable = nullptr
protectedinherited

◆ fUseNormalCorrection

G4bool G4VIntersectionLocator::fUseNormalCorrection = false
protectedinherited

◆ fVerboseLevel

G4int G4VIntersectionLocator::fVerboseLevel = 0
protectedinherited

◆ fWarnSteps

unsigned int G4MultiLevelLocator::fWarnSteps = 1000
private

Definition at line 84 of file G4MultiLevelLocator.hh.

Referenced by EstimateIntersectionPoint(), and SetWarnSteps().

◆ kCarTolerance

G4double G4VIntersectionLocator::kCarTolerance
protectedinherited

◆ max_depth

const G4int G4MultiLevelLocator::max_depth = 10
staticprivate

◆ ptrInterMedFT

G4FieldTrack* G4MultiLevelLocator::ptrInterMedFT[max_depth+1]
private

The documentation for this class was generated from the following files: