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