#include <G4SimpleLocator.hh>
Inheritance diagram for G4SimpleLocator:
Public Member Functions | |
G4SimpleLocator (G4Navigator *theNavigator) | |
~G4SimpleLocator () | |
G4bool | EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin) |
Definition at line 52 of file G4SimpleLocator.hh.
G4SimpleLocator::G4SimpleLocator | ( | G4Navigator * | theNavigator | ) |
Definition at line 39 of file G4SimpleLocator.cc.
00040 : G4VIntersectionLocator(theNavigator) 00041 { 00042 }
G4SimpleLocator::~G4SimpleLocator | ( | ) |
G4bool G4SimpleLocator::EstimateIntersectionPoint | ( | const G4FieldTrack & | curveStartPointTangent, | |
const G4FieldTrack & | curveEndPointTangent, | |||
const G4ThreeVector & | trialPoint, | |||
G4FieldTrack & | intersectPointTangent, | |||
G4bool & | recalculatedEndPoint, | |||
G4double & | fPreviousSafety, | |||
G4ThreeVector & | fPreviousSftOrigin | |||
) | [virtual] |
Implements G4VIntersectionLocator.
Definition at line 80 of file G4SimpleLocator.cc.
References G4VIntersectionLocator::AdjustmentOfFoundIntersection(), G4ChordFinder::ApproxCurvePointV(), FatalException, G4VIntersectionLocator::fiDeltaIntersection, G4VIntersectionLocator::fVerboseLevel, G4cout, G4endl, G4Exception(), G4VIntersectionLocator::GetAdjustementOfFoundIntersection(), G4VIntersectionLocator::GetChordFinderFor(), G4FieldTrack::GetCurveLength(), G4VIntersectionLocator::GetEpsilonStepFor(), G4FieldTrack::GetMomentumDir(), G4FieldTrack::GetMomentumDirection(), G4VIntersectionLocator::GetNavigatorFor(), G4FieldTrack::GetPosition(), G4VIntersectionLocator::GetSurfaceNormal(), G4VIntersectionLocator::IntersectChord(), JustWarning, G4VIntersectionLocator::kCarTolerance, G4Navigator::LocateGlobalPointWithinVolume(), G4VIntersectionLocator::printStatus(), G4VIntersectionLocator::ReEstimateEndpoint(), G4VIntersectionLocator::ReportTrialStep(), G4FieldTrack::SetPosition(), and sqr().
00088 { 00089 // Find Intersection Point ( A, B, E ) of true path AB - start at E. 00090 00091 G4bool found_approximate_intersection = false; 00092 G4bool there_is_no_intersection = false; 00093 00094 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity; 00095 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity; 00096 G4ThreeVector CurrentE_Point = TrialPoint; 00097 G4bool validNormalAtE = false; 00098 G4ThreeVector NormalAtEntry; 00099 00100 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct 00101 G4double NewSafety = 0.0; 00102 G4bool last_AF_intersection = false; 00103 G4bool final_section = true; // Shows whether current section is last 00104 // (i.e. B=full end) 00105 recalculatedEndPoint = false; 00106 00107 G4bool restoredFullEndpoint = false; 00108 00109 G4int substep_no = 0; 00110 00111 // Limits for substep number 00112 // 00113 const G4int max_substeps = 100000000; // Test 120 (old value 100 ) 00114 const G4int warn_substeps = 1000; // 100 00115 00116 // Statistics for substeps 00117 // 00118 static G4int max_no_seen= -1; 00119 00120 NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE); 00121 00122 #ifdef G4DEBUG_FIELD 00123 static G4double tolerance = 1.0e-8; 00124 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition(); 00125 if( (TrialPoint - StartPosition).mag() < tolerance * mm ) 00126 { 00127 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 00128 "GeomNav1002", JustWarning, 00129 "Intersection point F is exactly at start point A." ); 00130 } 00131 #endif 00132 00133 do 00134 { 00135 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 00136 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition(); 00137 00138 // F = a point on true AB path close to point E 00139 // (the closest if possible) 00140 // 00141 ApproxIntersecPointV = GetChordFinderFor() 00142 ->ApproxCurvePointV( CurrentA_PointVelocity, 00143 CurrentB_PointVelocity, 00144 CurrentE_Point, 00145 GetEpsilonStepFor()); 00146 // The above method is the key & most intuitive part ... 00147 00148 #ifdef G4DEBUG_FIELD 00149 if( ApproxIntersecPointV.GetCurveLength() > 00150 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) ) 00151 { 00152 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 00153 "GeomNav0003", FatalException, 00154 "Intermediate F point is past end B point!" ); 00155 } 00156 #endif 00157 00158 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition(); 00159 00160 // First check whether EF is small - then F is a good approx. point 00161 // Calculate the length and direction of the chord AF 00162 // 00163 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point; 00164 00165 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir(); 00166 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ; 00167 00168 G4ThreeVector ChordAB = Point_B - Point_A; 00169 00170 #ifdef G4DEBUG_FIELD 00171 G4VIntersectionLocator:: 00172 ReportTrialStep( substep_no, ChordAB, ChordEF_Vector, 00173 NewMomentumDir, NormalAtEntry, validNormalAtE ); 00174 #endif 00175 // Check Sign is always exiting !! TODO 00176 // Could ( > -epsilon) be used instead? 00177 // 00178 G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 ) 00179 || (! validNormalAtE ); // Invalid 00180 G4double EF_dist2= ChordEF_Vector.mag2(); 00181 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) ) 00182 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) ) 00183 { 00184 found_approximate_intersection = true; 00185 00186 // Create the "point" return value 00187 // 00188 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 00189 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point ); 00190 00191 if ( GetAdjustementOfFoundIntersection() ) 00192 { 00193 // Try to Get Correction of IntersectionPoint using SurfaceNormal() 00194 // 00195 G4ThreeVector IP; 00196 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection(); 00197 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A, 00198 CurrentE_Point, CurrentF_Point, MomentumDir, 00199 last_AF_intersection, IP, NewSafety, 00200 fPreviousSafety, fPreviousSftOrigin ); 00201 00202 if(goodCorrection) 00203 { 00204 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 00205 IntersectedOrRecalculatedFT.SetPosition(IP); 00206 } 00207 } 00208 00209 // Note: in order to return a point on the boundary, 00210 // we must return E. But it is F on the curve. 00211 // So we must "cheat": we are using the position at point E 00212 // and the velocity at point F !!! 00213 // 00214 // This must limit the length we can allow for displacement! 00215 } 00216 else // E is NOT close enough to the curve (ie point F) 00217 { 00218 // Check whether any volumes are encountered by the chord AF 00219 // --------------------------------------------------------- 00220 // First relocate to restore any Voxel etc information 00221 // in the Navigator before calling ComputeStep() 00222 // 00223 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A ); 00224 00225 G4ThreeVector PointG; // Candidate intersection point 00226 G4double stepLengthAF; 00227 G4bool usedNavigatorAF = false; 00228 G4bool Intersects_AF = IntersectChord( Point_A, 00229 CurrentF_Point, 00230 NewSafety, 00231 fPreviousSafety, 00232 fPreviousSftOrigin, 00233 stepLengthAF, 00234 PointG, 00235 &usedNavigatorAF ); 00236 last_AF_intersection = Intersects_AF; 00237 if( Intersects_AF ) 00238 { 00239 // G is our new Candidate for the intersection point. 00240 // It replaces "E" and we will repeat the test to see if 00241 // it is a good enough approximate point for us. 00242 // B <- F 00243 // E <- G 00244 00245 CurrentB_PointVelocity = ApproxIntersecPointV; 00246 CurrentE_Point = PointG; 00247 00248 // Need to recalculate the Exit Normal at the new PointG 00249 // Relies on a call to Navigator::ComputeStep in IntersectChord above 00250 // If the safety was adequate (for the step) this would NOT be called! 00251 // But this must not occur, no intersection can be found in that case, 00252 // so this branch, ie if( Intersects_AF ) would not be reached! 00253 // 00254 G4bool validNormalLast; 00255 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast ); 00256 validNormalAtE = validNormalLast; 00257 00258 // By moving point B, must take care if current 00259 // AF has no intersection to try current FB!! 00260 // 00261 final_section= false; 00262 00263 #ifdef G4VERBOSE 00264 if( fVerboseLevel > 3 ) 00265 { 00266 G4cout << "G4PiF::LI> Investigating intermediate point" 00267 << " at s=" << ApproxIntersecPointV.GetCurveLength() 00268 << " on way to full s=" 00269 << CurveEndPointVelocity.GetCurveLength() << G4endl; 00270 } 00271 #endif 00272 } 00273 else // not Intersects_AF 00274 { 00275 // In this case: 00276 // There is NO intersection of AF with a volume boundary. 00277 // We must continue the search in the segment FB! 00278 // 00279 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point ); 00280 00281 G4double stepLengthFB; 00282 G4ThreeVector PointH; 00283 G4bool usedNavigatorFB=false; 00284 00285 // Check whether any volumes are encountered by the chord FB 00286 // --------------------------------------------------------- 00287 00288 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 00289 NewSafety,fPreviousSafety, 00290 fPreviousSftOrigin, 00291 stepLengthFB, 00292 PointH, &usedNavigatorFB ); 00293 if( Intersects_FB ) 00294 { 00295 // There is an intersection of FB with a volume boundary 00296 // H <- First Intersection of Chord FB 00297 00298 // H is our new Candidate for the intersection point. 00299 // It replaces "E" and we will repeat the test to see if 00300 // it is a good enough approximate point for us. 00301 00302 // Note that F must be in volume volA (the same as A) 00303 // (otherwise AF would meet a volume boundary!) 00304 // A <- F 00305 // E <- H 00306 // 00307 CurrentA_PointVelocity = ApproxIntersecPointV; 00308 CurrentE_Point = PointH; 00309 00310 // Need to recalculate the Exit Normal at the new PointG 00311 // Relies on call to Navigator::ComputeStep in IntersectChord above 00312 // If safety was adequate (for the step) this would NOT be called! 00313 // But this must not occur, no intersection found in that case, 00314 // so this branch, i.e. if( Intersects_AF ) would not be reached! 00315 // 00316 G4bool validNormalLast; 00317 NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast ); 00318 validNormalAtE = validNormalLast; 00319 } 00320 else // not Intersects_FB 00321 { 00322 // There is NO intersection of FB with a volume boundary 00323 00324 if( final_section ) 00325 { 00326 // If B is the original endpoint, this means that whatever 00327 // volume(s) intersected the original chord, none touch the 00328 // smaller chords we have used. 00329 // The value of 'IntersectedOrRecalculatedFT' returned is 00330 // likely not valid 00331 00332 there_is_no_intersection = true; // real final_section 00333 } 00334 else 00335 { 00336 // We must restore the original endpoint 00337 00338 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B 00339 CurrentB_PointVelocity = CurveEndPointVelocity; 00340 restoredFullEndpoint = true; 00341 } 00342 } // Endif (Intersects_FB) 00343 } // Endif (Intersects_AF) 00344 00345 // Ensure that the new endpoints are not further apart in space 00346 // than on the curve due to different errors in the integration 00347 // 00348 G4double linDistSq, curveDist; 00349 linDistSq = ( CurrentB_PointVelocity.GetPosition() 00350 - CurrentA_PointVelocity.GetPosition() ).mag2(); 00351 curveDist = CurrentB_PointVelocity.GetCurveLength() 00352 - CurrentA_PointVelocity.GetCurveLength(); 00353 00354 // Change this condition for very strict parameters of propagation 00355 // 00356 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq ) 00357 { 00358 // Re-integrate to obtain a new B 00359 // 00360 G4FieldTrack newEndPointFT = 00361 ReEstimateEndpoint( CurrentA_PointVelocity, 00362 CurrentB_PointVelocity, 00363 linDistSq, // to avoid recalculation 00364 curveDist ); 00365 G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 00366 CurrentB_PointVelocity = newEndPointFT; 00367 00368 if( (final_section)) // real final section 00369 { 00370 recalculatedEndPoint = true; 00371 IntersectedOrRecalculatedFT = newEndPointFT; 00372 // So that we can return it, if it is the endpoint! 00373 } 00374 } 00375 if( curveDist < 0.0 ) 00376 { 00377 fVerboseLevel = 5; // Print out a maximum of information 00378 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 00379 -1.0, NewSafety, substep_no ); 00380 std::ostringstream message; 00381 message << "Error in advancing propagation." << G4endl 00382 << " Point A (start) is " << CurrentA_PointVelocity 00383 << G4endl 00384 << " Point B (end) is " << CurrentB_PointVelocity 00385 << G4endl 00386 << " Curve distance is " << curveDist << G4endl 00387 << G4endl 00388 << "The final curve point is not further along" 00389 << " than the original!" << G4endl; 00390 00391 if( recalculatedEndPoint ) 00392 { 00393 message << "Recalculation of EndPoint was called with fEpsStep= " 00394 << GetEpsilonStepFor() << G4endl; 00395 } 00396 message.precision(20); 00397 message << " Point A (Curve start) is " << CurveStartPointVelocity 00398 << G4endl 00399 << " Point B (Curve end) is " << CurveEndPointVelocity 00400 << G4endl 00401 << " Point A (Current start) is " << CurrentA_PointVelocity 00402 << G4endl 00403 << " Point B (Current end) is " << CurrentB_PointVelocity 00404 << G4endl 00405 << " Point E (Trial Point) is " << CurrentE_Point 00406 << G4endl 00407 << " Point F (Intersection) is " << ApproxIntersecPointV 00408 << G4endl 00409 << " LocateIntersection parameters are : Substep no= " 00410 << substep_no; 00411 00412 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 00413 "GeomNav0003", FatalException, message); 00414 } 00415 00416 if(restoredFullEndpoint) 00417 { 00418 final_section = restoredFullEndpoint; 00419 restoredFullEndpoint = false; 00420 } 00421 } // EndIf ( E is close enough to the curve, ie point F. ) 00422 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 00423 00424 #ifdef G4DEBUG_LOCATE_INTERSECTION 00425 static G4int trigger_substepno_print= warn_substeps - 20; 00426 00427 if( substep_no >= trigger_substepno_print ) 00428 { 00429 G4cout << "Difficulty in converging in " 00430 << "G4SimpleLocator::EstimateIntersectionPoint():" 00431 << G4endl 00432 << " Substep no = " << substep_no << G4endl; 00433 if( substep_no == trigger_substepno_print ) 00434 { 00435 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 00436 -1.0, NewSafety, 0); 00437 } 00438 G4cout << " State of point A: "; 00439 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 00440 -1.0, NewSafety, substep_no-1, 0); 00441 G4cout << " State of point B: "; 00442 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 00443 -1.0, NewSafety, substep_no); 00444 } 00445 #endif 00446 substep_no++; 00447 00448 } while ( ( ! found_approximate_intersection ) 00449 && ( ! there_is_no_intersection ) 00450 && ( substep_no <= max_substeps) ); // UNTIL found or failed 00451 00452 if( substep_no > max_no_seen ) 00453 { 00454 max_no_seen = substep_no; 00455 #ifdef G4DEBUG_LOCATE_INTERSECTION 00456 if( max_no_seen > warn_substeps ) 00457 { 00458 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 00459 } 00460 #endif 00461 } 00462 00463 if( ( substep_no >= max_substeps) 00464 && !there_is_no_intersection 00465 && !found_approximate_intersection ) 00466 { 00467 G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl 00468 << " Start and Endpoint of Requested Step:" << G4endl; 00469 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 00470 -1.0, NewSafety, 0); 00471 G4cout << G4endl 00472 << " Start and end-point of current Sub-Step:" << G4endl; 00473 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 00474 -1.0, NewSafety, substep_no-1); 00475 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 00476 -1.0, NewSafety, substep_no); 00477 00478 std::ostringstream message; 00479 message << "Convergence is requiring too many substeps: " 00480 << substep_no << G4endl 00481 << " Abandoning effort to intersect." << G4endl 00482 << " Found intersection = " 00483 << found_approximate_intersection << G4endl 00484 << " Intersection exists = " 00485 << !there_is_no_intersection << G4endl; 00486 message.precision(10); 00487 G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 00488 G4double full_len = CurveEndPointVelocity.GetCurveLength(); 00489 message << " Undertaken only length: " << done_len 00490 << " out of " << full_len << " required." << G4endl 00491 << " Remaining length = " << full_len-done_len; 00492 00493 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 00494 "GeomNav0003", FatalException, message); 00495 } 00496 else if( substep_no >= warn_substeps ) 00497 { 00498 std::ostringstream message; 00499 message.precision(10); 00500 00501 message << "Many substeps while trying to locate intersection." << G4endl 00502 << " Undertaken length: " 00503 << CurrentB_PointVelocity.GetCurveLength() 00504 << " - Needed: " << substep_no << " substeps." << G4endl 00505 << " Warning level = " << warn_substeps 00506 << " and maximum substeps = " << max_substeps; 00507 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 00508 "GeomNav1002", JustWarning, message); 00509 } 00510 return !there_is_no_intersection; // Success or failure 00511 }