G4BrentLocator Class Reference

#include <G4BrentLocator.hh>

Inheritance diagram for G4BrentLocator:

G4VIntersectionLocator

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)

Detailed Description

Definition at line 50 of file G4BrentLocator.hh.


Constructor & Destructor Documentation

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 }


Member Function Documentation

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 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:33 2013 for Geant4 by  doxygen 1.4.7