#include <G4BrentLocator.hh>
Inheritance diagram for G4BrentLocator:
Public Member Functions | |
G4BrentLocator (G4Navigator *theNavigator) | |
~G4BrentLocator () | |
G4bool | EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin) |
Definition at line 50 of file G4BrentLocator.hh.
G4BrentLocator::G4BrentLocator | ( | G4Navigator * | theNavigator | ) |
Definition at line 38 of file G4BrentLocator.cc.
00039 : G4VIntersectionLocator(theNavigator) 00040 { 00041 // In case of too slow progress in finding Intersection Point 00042 // intermediates Points on the Track must be stored. 00043 // Initialise the array of Pointers [max_depth+1] to do this 00044 00045 G4ThreeVector zeroV(0.0,0.0,0.0); 00046 for (G4int idepth=0; idepth<max_depth+1; idepth++ ) 00047 { 00048 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.); 00049 } 00050 00051 // Counters for Locator 00052 00053 // Counter for Maximum Number Of Trial before Intersection Found 00054 // 00055 maxNumberOfStepsForIntersection=0; 00056 00057 // Counter for Number Of Calls to ReIntegrationEndPoint Method 00058 // 00059 maxNumberOfCallsToReIntegration=0; 00060 maxNumberOfCallsToReIntegration_depth=0; 00061 }
G4BrentLocator::~G4BrentLocator | ( | ) |
Definition at line 63 of file G4BrentLocator.cc.
References G4VIntersectionLocator::fVerboseLevel, G4cout, and G4endl.
00064 { 00065 for ( G4int idepth=0; idepth<max_depth+1; idepth++) 00066 { 00067 delete ptrInterMedFT[idepth]; 00068 } 00069 #ifdef G4DEBUG_FIELD 00070 if(fVerboseLevel>0) 00071 { 00072 G4cout << "G4BrentLocator::Location with Max Number of Steps=" 00073 << maxNumberOfStepsForIntersection<<G4endl; 00074 G4cout << "G4BrentLocator::ReIntegrateEndPoint was called " 00075 << maxNumberOfCallsToReIntegration 00076 << " times and for depth algorithm " 00077 << maxNumberOfCallsToReIntegration_depth << " times." << G4endl; 00078 } 00079 #endif 00080 }
G4bool G4BrentLocator::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 G4BrentLocator.cc.
References G4MagInt_Driver::AccurateAdvance(), G4VIntersectionLocator::AdjustmentOfFoundIntersection(), G4ChordFinder::ApproxCurvePointS(), G4ChordFinder::ApproxCurvePointV(), FatalException, G4VIntersectionLocator::fiDeltaIntersection, G4VIntersectionLocator::fVerboseLevel, G4cerr, G4cout, G4endl, G4Exception(), G4VIntersectionLocator::GetAdjustementOfFoundIntersection(), G4VIntersectionLocator::GetChordFinderFor(), G4FieldTrack::GetCurveLength(), G4VIntersectionLocator::GetDeltaIntersectionFor(), G4VIntersectionLocator::GetEpsilonStepFor(), G4ChordFinder::GetIntegrationDriver(), 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().
00124 { 00125 // Find Intersection Point ( A, B, E ) of true path AB - start at E. 00126 00127 G4bool found_approximate_intersection = false; 00128 G4bool there_is_no_intersection = false; 00129 00130 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity; 00131 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity; 00132 G4ThreeVector CurrentE_Point = TrialPoint; 00133 G4bool validNormalAtE = false; 00134 G4ThreeVector NormalAtEntry; 00135 00136 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct 00137 G4double NewSafety = 0.0; 00138 G4bool last_AF_intersection = false; 00139 00140 // G4bool final_section= true; // Shows whether current section is last 00141 // (i.e. B=full end) 00142 G4bool first_section = true; 00143 recalculatedEndPoint = false; 00144 00145 G4bool restoredFullEndpoint = false; 00146 00147 G4int oldprc; // cout, cerr precision 00148 G4int substep_no = 0; 00149 00150 // Limits for substep number 00151 // 00152 const G4int max_substeps= 10000; // Test 120 (old value 100 ) 00153 const G4int warn_substeps= 1000; // 100 00154 00155 // Statistics for substeps 00156 // 00157 static G4int max_no_seen= -1; 00158 00159 // Counter for restarting Bintermed 00160 // 00161 G4int restartB = 0; 00162 00163 //-------------------------------------------------------------------------- 00164 // Algorithm for the case if progress in founding intersection is too slow. 00165 // Process is defined too slow if after N=param_substeps advances on the 00166 // path, it will be only 'fraction_done' of the total length. 00167 // In this case the remaining length is divided in two half and 00168 // the loop is restarted for each half. 00169 // If progress is still too slow, the division in two halfs continue 00170 // until 'max_depth'. 00171 //-------------------------------------------------------------------------- 00172 00173 const G4int param_substeps=50; // Test value for the maximum number 00174 // of substeps 00175 const G4double fraction_done=0.3; 00176 00177 G4bool Second_half = false; // First half or second half of divided step 00178 00179 NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE); 00180 00181 // We need to know this for the 'final_section': 00182 // real 'final_section' or first half 'final_section' 00183 // In algorithm it is considered that the 'Second_half' is true 00184 // and it becomes false only if we are in the first-half of level 00185 // depthness or if we are in the first section 00186 00187 G4int depth=0; // Depth counts how many subdivisions of initial step made 00188 00189 #ifdef G4DEBUG_FIELD 00190 static G4double tolerance= 1.0e-8; 00191 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition(); 00192 if( (TrialPoint - StartPosition).mag() < tolerance * mm ) 00193 { 00194 G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 00195 "GeomNav1002", JustWarning, 00196 "Intersection point F is exactly at start point A." ); 00197 } 00198 #endif 00199 00200 // Intermediates Points on the Track = Subdivided Points must be stored. 00201 // Give the initial values to 'InterMedFt' 00202 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint' 00203 // 00204 *ptrInterMedFT[0] = CurveEndPointVelocity; 00205 for (G4int idepth=1; idepth<max_depth+1; idepth++ ) 00206 { 00207 *ptrInterMedFT[idepth]=CurveStartPointVelocity; 00208 } 00209 00210 //Final_section boolean store 00211 G4bool fin_section_depth[max_depth]; 00212 for (G4int idepth=0; idepth<max_depth; idepth++ ) 00213 { 00214 fin_section_depth[idepth]=true; 00215 } 00216 00217 // 'SubStartPoint' is needed to calculate the length of the divided step 00218 // 00219 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity; 00220 00221 do 00222 { 00223 G4int substep_no_p = 0; 00224 G4bool sub_final_section = false; // the same as final_section, 00225 // but for 'sub_section' 00226 SubStart_PointVelocity = CurrentA_PointVelocity; 00227 do // REPEAT param 00228 { 00229 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 00230 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition(); 00231 00232 // F = a point on true AB path close to point E 00233 // (the closest if possible) 00234 // 00235 if(substep_no_p==0) 00236 { 00237 ApproxIntersecPointV = GetChordFinderFor() 00238 ->ApproxCurvePointV( CurrentA_PointVelocity, 00239 CurrentB_PointVelocity, 00240 CurrentE_Point, 00241 GetEpsilonStepFor()); 00242 // The above method is the key & most intuitive part ... 00243 } 00244 #ifdef G4DEBUG_FIELD 00245 if( ApproxIntersecPointV.GetCurveLength() > 00246 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) ) 00247 { 00248 G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 00249 "GeomNav0003", FatalException, 00250 "Intermediate F point is past end B point" ); 00251 } 00252 #endif 00253 00254 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition(); 00255 00256 // First check whether EF is small - then F is a good approx. point 00257 // Calculate the length and direction of the chord AF 00258 // 00259 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point; 00260 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir(); 00261 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ; 00262 00263 #ifdef G4DEBUG_FIELD 00264 G4ThreeVector ChordAB = Point_B - Point_A; 00265 00266 G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB, 00267 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE ); 00268 #endif 00269 00270 G4bool adequate_angle; 00271 adequate_angle = ( MomDir_dot_Norm >= 0.0 ) // Can use -epsilon instead. 00272 || (! validNormalAtE ); // Makes criterion invalid 00273 G4double EF_dist2 = ChordEF_Vector.mag2(); 00274 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) ) 00275 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) ) 00276 { 00277 found_approximate_intersection = true; 00278 00279 // Create the "point" return value 00280 // 00281 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 00282 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point ); 00283 00284 if ( GetAdjustementOfFoundIntersection() ) 00285 { 00286 // Try to Get Correction of IntersectionPoint using SurfaceNormal() 00287 // 00288 G4ThreeVector IP; 00289 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection(); 00290 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A, 00291 CurrentE_Point, CurrentF_Point, MomentumDir, 00292 last_AF_intersection, IP, NewSafety, 00293 fPreviousSafety, fPreviousSftOrigin ); 00294 if ( goodCorrection ) 00295 { 00296 IntersectedOrRecalculatedFT = ApproxIntersecPointV; 00297 IntersectedOrRecalculatedFT.SetPosition(IP); 00298 } 00299 } 00300 00301 // Note: in order to return a point on the boundary, 00302 // we must return E. But it is F on the curve. 00303 // So we must "cheat": we are using the position at point E 00304 // and the velocity at point F !!! 00305 // 00306 // This must limit the length we can allow for displacement! 00307 } 00308 else // E is NOT close enough to the curve (ie point F) 00309 { 00310 // Check whether any volumes are encountered by the chord AF 00311 // --------------------------------------------------------- 00312 // First relocate to restore any Voxel etc information 00313 // in the Navigator before calling ComputeStep() 00314 // 00315 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A ); 00316 00317 G4ThreeVector PointG; // Candidate intersection point 00318 G4double stepLengthAF; 00319 G4bool usedNavigatorAF = false; 00320 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point, 00321 NewSafety,fPreviousSafety, 00322 fPreviousSftOrigin, 00323 stepLengthAF, 00324 PointG, 00325 &usedNavigatorAF); 00326 last_AF_intersection = Intersects_AF; 00327 if( Intersects_AF ) 00328 { 00329 // G is our new Candidate for the intersection point. 00330 // It replaces "E" and we will repeat the test to see if 00331 // it is a good enough approximate point for us. 00332 // B <- F 00333 // E <- G 00334 // 00335 G4FieldTrack EndPoint = ApproxIntersecPointV; 00336 ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS( 00337 CurrentA_PointVelocity, CurrentB_PointVelocity, 00338 EndPoint,CurrentE_Point, CurrentF_Point,PointG, 00339 true, GetEpsilonStepFor() ); 00340 CurrentB_PointVelocity = EndPoint; 00341 CurrentE_Point = PointG; 00342 00343 // Need to recalculate the Exit Normal at the new PointG 00344 // Know that a call was made to Navigator::ComputeStep in 00345 // IntersectChord above. 00346 // 00347 G4bool validNormalLast; 00348 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast ); 00349 validNormalAtE = validNormalLast; 00350 00351 // By moving point B, must take care if current 00352 // AF has no intersection to try current FB!! 00353 // 00354 fin_section_depth[depth] = false; 00355 #ifdef G4VERBOSE 00356 if( fVerboseLevel > 3 ) 00357 { 00358 G4cout << "G4PiF::LI> Investigating intermediate point" 00359 << " at s=" << ApproxIntersecPointV.GetCurveLength() 00360 << " on way to full s=" 00361 << CurveEndPointVelocity.GetCurveLength() << G4endl; 00362 } 00363 #endif 00364 } 00365 else // not Intersects_AF 00366 { 00367 // In this case: 00368 // There is NO intersection of AF with a volume boundary. 00369 // We must continue the search in the segment FB! 00370 // 00371 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point ); 00372 00373 G4double stepLengthFB; 00374 G4ThreeVector PointH; 00375 G4bool usedNavigatorFB = false; 00376 00377 // Check whether any volumes are encountered by the chord FB 00378 // --------------------------------------------------------- 00379 00380 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 00381 NewSafety,fPreviousSafety, 00382 fPreviousSftOrigin, 00383 stepLengthFB, 00384 PointH, 00385 &usedNavigatorFB); 00386 if( Intersects_FB ) 00387 { 00388 // There is an intersection of FB with a volume boundary 00389 // H <- First Intersection of Chord FB 00390 00391 // H is our new Candidate for the intersection point. 00392 // It replaces "E" and we will repeat the test to see if 00393 // it is a good enough approximate point for us. 00394 00395 // Note that F must be in volume volA (the same as A) 00396 // (otherwise AF would meet a volume boundary!) 00397 // A <- F 00398 // E <- H 00399 // 00400 G4FieldTrack InterMed=ApproxIntersecPointV; 00401 ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS( 00402 CurrentA_PointVelocity,CurrentB_PointVelocity, 00403 InterMed,CurrentE_Point,CurrentF_Point,PointH, 00404 false,GetEpsilonStepFor()); 00405 CurrentA_PointVelocity = InterMed; 00406 CurrentE_Point = PointH; 00407 00408 // Need to recalculate the Exit Normal at the new PointG 00409 // 00410 G4bool validNormalLast; 00411 NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast ); 00412 validNormalAtE= validNormalLast; 00413 } 00414 else // not Intersects_FB 00415 { 00416 // There is NO intersection of FB with a volume boundary 00417 00418 if( fin_section_depth[depth] ) 00419 { 00420 // If B is the original endpoint, this means that whatever 00421 // volume(s) intersected the original chord, none touch the 00422 // smaller chords we have used. 00423 // The value of 'IntersectedOrRecalculatedFT' returned is 00424 // likely not valid 00425 00426 // Check on real final_section or SubEndSection 00427 // 00428 if( ((Second_half)&&(depth==0)) || (first_section) ) 00429 { 00430 there_is_no_intersection = true; // real final_section 00431 } 00432 else 00433 { 00434 // end of subsection, not real final section 00435 // exit from the and go to the depth-1 level 00436 00437 substep_no_p = param_substeps+2; // exit from the loop 00438 00439 // but 'Second_half' is still true because we need to find 00440 // the 'CurrentE_point' for the next loop 00441 // 00442 Second_half = true; 00443 sub_final_section = true; 00444 } 00445 } 00446 else 00447 { 00448 if(depth==0) 00449 { 00450 // We must restore the original endpoint 00451 // 00452 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B 00453 CurrentB_PointVelocity = CurveEndPointVelocity; 00454 SubStart_PointVelocity = CurrentA_PointVelocity; 00455 ApproxIntersecPointV = GetChordFinderFor() 00456 ->ApproxCurvePointV( CurrentA_PointVelocity, 00457 CurrentB_PointVelocity, 00458 CurrentE_Point, 00459 GetEpsilonStepFor()); 00460 00461 restoredFullEndpoint = true; 00462 restartB++; // counter 00463 } 00464 else 00465 { 00466 // We must restore the depth endpoint 00467 // 00468 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B 00469 CurrentB_PointVelocity = *ptrInterMedFT[depth]; 00470 SubStart_PointVelocity = CurrentA_PointVelocity; 00471 ApproxIntersecPointV = GetChordFinderFor() 00472 ->ApproxCurvePointV( CurrentA_PointVelocity, 00473 CurrentB_PointVelocity, 00474 CurrentE_Point, 00475 GetEpsilonStepFor()); 00476 restoredFullEndpoint = true; 00477 restartB++; // counter 00478 } 00479 } 00480 } // Endif (Intersects_FB) 00481 } // Endif (Intersects_AF) 00482 00483 // Ensure that the new endpoints are not further apart in space 00484 // than on the curve due to different errors in the integration 00485 // 00486 G4double linDistSq, curveDist; 00487 linDistSq = ( CurrentB_PointVelocity.GetPosition() 00488 - CurrentA_PointVelocity.GetPosition() ).mag2(); 00489 curveDist = CurrentB_PointVelocity.GetCurveLength() 00490 - CurrentA_PointVelocity.GetCurveLength(); 00491 00492 // Change this condition for very strict parameters of propagation 00493 // 00494 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq ) 00495 { 00496 // Re-integrate to obtain a new B 00497 // 00498 G4FieldTrack newEndPointFT= 00499 ReEstimateEndpoint( CurrentA_PointVelocity, 00500 CurrentB_PointVelocity, 00501 linDistSq, // to avoid recalculation 00502 curveDist ); 00503 G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 00504 CurrentB_PointVelocity = newEndPointFT; 00505 00506 if ( (fin_section_depth[depth]) // real final section 00507 &&( first_section || ((Second_half)&&(depth==0)) ) ) 00508 { 00509 recalculatedEndPoint = true; 00510 IntersectedOrRecalculatedFT = newEndPointFT; 00511 // So that we can return it, if it is the endpoint! 00512 } 00513 } 00514 if( curveDist < 0.0 ) 00515 { 00516 fVerboseLevel = 5; // Print out a maximum of information 00517 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 00518 -1.0, NewSafety, substep_no ); 00519 std::ostringstream message; 00520 message << "Error in advancing propagation." << G4endl 00521 << " Error in advancing propagation." << G4endl 00522 << " Point A (start) is " << CurrentA_PointVelocity 00523 << G4endl 00524 << " Point B (end) is " << CurrentB_PointVelocity 00525 << G4endl 00526 << " Curve distance is " << curveDist << G4endl 00527 << G4endl 00528 << "The final curve point is not further along" 00529 << " than the original!" << G4endl; 00530 00531 if( recalculatedEndPoint ) 00532 { 00533 message << "Recalculation of EndPoint was called with fEpsStep= " 00534 << GetEpsilonStepFor() << G4endl; 00535 } 00536 oldprc = G4cerr.precision(20); 00537 message << " Point A (Curve start) is " << CurveStartPointVelocity 00538 << G4endl 00539 << " Point B (Curve end) is " << CurveEndPointVelocity 00540 << G4endl 00541 << " Point A (Current start) is " << CurrentA_PointVelocity 00542 << G4endl 00543 << " Point B (Current end) is " << CurrentB_PointVelocity 00544 << G4endl 00545 << " Point S (Sub start) is " << SubStart_PointVelocity 00546 << G4endl 00547 << " Point E (Trial Point) is " << CurrentE_Point 00548 << G4endl 00549 << " Old Point F(Intersection) is " << CurrentF_Point 00550 << G4endl 00551 << " New Point F(Intersection) is " << ApproxIntersecPointV 00552 << G4endl 00553 << " LocateIntersection parameters are : Substep no= " 00554 << substep_no << G4endl 00555 << " Substep depth no= "<< substep_no_p << " Depth= " 00556 << depth << G4endl 00557 << " Restarted no= "<< restartB << " Epsilon= " 00558 << GetEpsilonStepFor() <<" DeltaInters= " 00559 << GetDeltaIntersectionFor(); 00560 G4cerr.precision( oldprc ); 00561 00562 G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 00563 "GeomNav0003", FatalException, message); 00564 } 00565 00566 if(restoredFullEndpoint) 00567 { 00568 fin_section_depth[depth] = restoredFullEndpoint; 00569 restoredFullEndpoint = false; 00570 } 00571 } // EndIf ( E is close enough to the curve, ie point F. ) 00572 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 00573 00574 #ifdef G4DEBUG_LOCATE_INTERSECTION 00575 static G4int trigger_substepno_print= warn_substeps - 20 ; 00576 00577 if( substep_no >= trigger_substepno_print ) 00578 { 00579 G4cout << "Difficulty in converging in " 00580 << "G4BrentLocator::EstimateIntersectionPoint()" 00581 << G4endl 00582 << " Substep no = " << substep_no << G4endl; 00583 if( substep_no == trigger_substepno_print ) 00584 { 00585 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 00586 -1.0, NewSafety, 0); 00587 } 00588 G4cout << " State of point A: "; 00589 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 00590 -1.0, NewSafety, substep_no-1, 0); 00591 G4cout << " State of point B: "; 00592 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 00593 -1.0, NewSafety, substep_no); 00594 } 00595 #endif 00596 substep_no++; 00597 substep_no_p++; 00598 00599 } while ( ( ! found_approximate_intersection ) 00600 && ( ! there_is_no_intersection ) 00601 && ( substep_no_p <= param_substeps) ); // UNTIL found or 00602 // failed param substep 00603 first_section = false; 00604 00605 if( (!found_approximate_intersection) && (!there_is_no_intersection) ) 00606 { 00607 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength() 00608 - SubStart_PointVelocity.GetCurveLength()); 00609 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength() 00610 - SubStart_PointVelocity.GetCurveLength()); 00611 00612 G4double stepLengthAB; 00613 G4ThreeVector PointGe; 00614 00615 // Check if progress is too slow and if it possible to go deeper, 00616 // then halve the step if so 00617 // 00618 if ( ( did_len < fraction_done*all_len ) 00619 && (depth<max_depth) && (!sub_final_section) ) 00620 { 00621 Second_half=false; 00622 depth++; 00623 00624 G4double Sub_len = (all_len-did_len)/(2.); 00625 G4FieldTrack start = CurrentA_PointVelocity; 00626 G4MagInt_Driver* integrDriver = 00627 GetChordFinderFor()->GetIntegrationDriver(); 00628 integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor()); 00629 *ptrInterMedFT[depth] = start; 00630 CurrentB_PointVelocity = *ptrInterMedFT[depth]; 00631 00632 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop 00633 // 00634 SubStart_PointVelocity = CurrentA_PointVelocity; 00635 00636 // Find new trial intersection point needed at start of the loop 00637 // 00638 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 00639 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition(); 00640 00641 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A); 00642 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, 00643 NewSafety, fPreviousSafety, 00644 fPreviousSftOrigin,stepLengthAB, 00645 PointGe); 00646 if( Intersects_AB ) 00647 { 00648 last_AF_intersection = Intersects_AB; 00649 CurrentE_Point = PointGe; 00650 fin_section_depth[depth]=true; 00651 00652 // Need to recalculate the Exit Normal at the new PointG 00653 // 00654 G4bool validNormalAB; 00655 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB ); 00656 validNormalAtE= validNormalAB; 00657 } 00658 else 00659 { 00660 // No intersection found for first part of curve 00661 // (CurrentA,InterMedPoint[depth]). Go to the second part 00662 // 00663 Second_half = true; 00664 } 00665 } // if did_len 00666 00667 if( (Second_half)&&(depth!=0) ) 00668 { 00669 // Second part of curve (InterMed[depth],Intermed[depth-1]) ) 00670 // On the depth-1 level normally we are on the 'second_half' 00671 00672 Second_half = true; 00673 00674 // Find new trial intersection point needed at start of the loop 00675 // 00676 SubStart_PointVelocity = *ptrInterMedFT[depth]; 00677 CurrentA_PointVelocity = *ptrInterMedFT[depth]; 00678 CurrentB_PointVelocity = *ptrInterMedFT[depth-1]; 00679 // Ensure that the new endpoints are not further apart in space 00680 // than on the curve due to different errors in the integration 00681 // 00682 G4double linDistSq, curveDist; 00683 linDistSq = ( CurrentB_PointVelocity.GetPosition() 00684 - CurrentA_PointVelocity.GetPosition() ).mag2(); 00685 curveDist = CurrentB_PointVelocity.GetCurveLength() 00686 - CurrentA_PointVelocity.GetCurveLength(); 00687 if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq ) 00688 { 00689 // Re-integrate to obtain a new B 00690 // 00691 G4FieldTrack newEndPointFT= 00692 ReEstimateEndpoint( CurrentA_PointVelocity, 00693 CurrentB_PointVelocity, 00694 linDistSq, // to avoid recalculation 00695 curveDist ); 00696 G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 00697 CurrentB_PointVelocity = newEndPointFT; 00698 if (depth==1) 00699 { 00700 recalculatedEndPoint = true; 00701 IntersectedOrRecalculatedFT = newEndPointFT; 00702 // So that we can return it, if it is the endpoint! 00703 } 00704 } 00705 00706 00707 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 00708 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition(); 00709 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A); 00710 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety, 00711 fPreviousSafety, 00712 fPreviousSftOrigin,stepLengthAB, PointGe); 00713 if( Intersects_AB ) 00714 { 00715 last_AF_intersection = Intersects_AB; 00716 CurrentE_Point = PointGe; 00717 00718 G4bool validNormalAB; 00719 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB ); 00720 validNormalAtE = validNormalAB; 00721 } 00722 00723 depth--; 00724 fin_section_depth[depth]=true; 00725 } 00726 } // if(!found_aproximate_intersection) 00727 00728 } while ( ( ! found_approximate_intersection ) 00729 && ( ! there_is_no_intersection ) 00730 && ( substep_no <= max_substeps) ); // UNTIL found or failed 00731 00732 if( substep_no > max_no_seen ) 00733 { 00734 max_no_seen = substep_no; 00735 #ifdef G4DEBUG_LOCATE_INTERSECTION 00736 if( max_no_seen > warn_substeps ) 00737 { 00738 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 00739 } 00740 #endif 00741 } 00742 00743 if( ( substep_no >= max_substeps) 00744 && !there_is_no_intersection 00745 && !found_approximate_intersection ) 00746 { 00747 G4cout << "ERROR - G4BrentLocator::EstimateIntersectionPoint()" << G4endl 00748 << " Start and end-point of requested Step:" << G4endl; 00749 printStatus( CurveStartPointVelocity, CurveEndPointVelocity, 00750 -1.0, NewSafety, 0); 00751 G4cout << " Start and end-point of current Sub-Step:" << G4endl; 00752 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity, 00753 -1.0, NewSafety, substep_no-1); 00754 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity, 00755 -1.0, NewSafety, substep_no); 00756 std::ostringstream message; 00757 message << "Too many substeps!" << G4endl 00758 << " Convergence is requiring too many substeps: " 00759 << substep_no << G4endl 00760 << " Abandoning effort to intersect. " << G4endl 00761 << " Found intersection = " 00762 << found_approximate_intersection << G4endl 00763 << " Intersection exists = " 00764 << !there_is_no_intersection << G4endl; 00765 oldprc = G4cout.precision( 10 ); 00766 G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 00767 G4double full_len = CurveEndPointVelocity.GetCurveLength(); 00768 message << " Undertaken only length: " << done_len 00769 << " out of " << full_len << " required." << G4endl 00770 << " Remaining length = " << full_len - done_len; 00771 G4cout.precision( oldprc ); 00772 00773 G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 00774 "GeomNav0003", FatalException, message); 00775 } 00776 else if( substep_no >= warn_substeps ) 00777 { 00778 oldprc= G4cout.precision( 10 ); 00779 std::ostringstream message; 00780 message << "Many substeps while trying to locate intersection." 00781 << G4endl 00782 << " Undertaken length: " 00783 << CurrentB_PointVelocity.GetCurveLength() 00784 << " - Needed: " << substep_no << " substeps." << G4endl 00785 << " Warning level = " << warn_substeps 00786 << " and maximum substeps = " << max_substeps; 00787 G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 00788 "GeomNav1002", JustWarning, message); 00789 G4cout.precision( oldprc ); 00790 } 00791 return !there_is_no_intersection; // Success or failure 00792 }