G4MultiLevelLocator.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Class G4MultiLevelLocator implementation
00029 //
00030 // 27.10.08 - Tatiana Nikitina.
00031 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
00032 // ---------------------------------------------------------------------------
00033 
00034 #include <iomanip>
00035 
00036 #include "G4ios.hh"
00037 #include "G4MultiLevelLocator.hh"
00038 
00039 G4MultiLevelLocator::G4MultiLevelLocator(G4Navigator *theNavigator)
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 }      
00052 
00053 G4MultiLevelLocator::~G4MultiLevelLocator()
00054 {
00055   for ( G4int idepth=0; idepth<max_depth+1; idepth++)
00056   {
00057     delete ptrInterMedFT[idepth];
00058   }
00059 }
00060 
00061 // --------------------------------------------------------------------------
00062 // G4bool G4PropagatorInField::LocateIntersectionPoint( 
00063 //        const G4FieldTrack&       CurveStartPointVelocity,   // A
00064 //        const G4FieldTrack&       CurveEndPointVelocity,     // B
00065 //        const G4ThreeVector&      TrialPoint,                // E
00066 //              G4FieldTrack&       IntersectedOrRecalculated  // Output
00067 //              G4bool&             recalculated )             // Out
00068 // --------------------------------------------------------------------------
00069 //
00070 // Function that returns the intersection of the true path with the surface
00071 // of the current volume (either the external one or the inner one with one
00072 // of the daughters:
00073 //
00074 //     A = Initial point
00075 //     B = another point 
00076 //
00077 // Both A and B are assumed to be on the true path:
00078 //
00079 //     E is the first point of intersection of the chord AB with 
00080 //     a volume other than A  (on the surface of A or of a daughter)
00081 //
00082 // Convention of Use :
00083 //     i) If it returns "true", then IntersectionPointVelocity is set
00084 //       to the approximate intersection point.
00085 //    ii) If it returns "false", no intersection was found.
00086 //          The validity of IntersectedOrRecalculated depends on 'recalculated'
00087 //        a) if latter is false, then IntersectedOrRecalculated is invalid. 
00088 //        b) if latter is true,  then IntersectedOrRecalculated is
00089 //             the new endpoint, due to a re-integration.
00090 // --------------------------------------------------------------------------
00091 // NOTE: implementation taken from G4PropagatorInField
00092 //
00093 G4bool G4MultiLevelLocator::EstimateIntersectionPoint( 
00094          const  G4FieldTrack&       CurveStartPointVelocity,       // A
00095          const  G4FieldTrack&       CurveEndPointVelocity,         // B
00096          const  G4ThreeVector&      TrialPoint,                    // E
00097                 G4FieldTrack&       IntersectedOrRecalculatedFT,   // Output
00098                 G4bool&             recalculatedEndPoint,          // Out
00099                 G4double&           previousSafety,                // In/Out
00100                 G4ThreeVector&      previousSftOrigin)             // In/Out
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 }

Generated on Mon May 27 17:48:54 2013 for Geant4 by  doxygen 1.4.7