Ignore:
Timestamp:
Feb 16, 2009, 10:14:30 AM (15 years ago)
Author:
garnier
Message:

en test de gl2ps. Problemes de libraries

Location:
trunk/source/geometry/navigation/src
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/navigation/src/G4AuxiliaryNavServices.cc

    r850 r921  
    2525//
    2626// $Id: G4AuxiliaryNavServices.cc,v 1.3 2006/06/29 18:36:32 gunter Exp $
    27 // GEANT4 tag $Name: HEAD $
     27// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2828//
    2929// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4DrawVoxels.cc

    r850 r921  
    2626//
    2727// $Id: G4DrawVoxels.cc,v 1.4 2006/06/29 18:36:34 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/src/G4ErrorPropagationNavigator.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4ErrorPropagationNavigator.cc,v 1.1 2007/05/16 12:49:18 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4ErrorPropagationNavigator.cc,v 1.2 2008/10/24 14:00:03 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    126126G4double G4ErrorPropagationNavigator::
    127127ComputeSafety( const G4ThreeVector &pGlobalpoint,
    128                const G4double pMaxLength )
     128               const G4double pMaxLength,
     129               const G4bool keepState )
    129130{
    130   G4double newSafety = G4Navigator::ComputeSafety(pGlobalpoint, pMaxLength);
     131  G4double newSafety = G4Navigator::ComputeSafety(pGlobalpoint,
     132                                                  pMaxLength, keepState);
    131133
    132134  G4ErrorPropagatorData *g4edata

  • trunk/source/geometry/navigation/src/G4GeomTestErrorList.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestErrorList.cc,v 1.3 2006/06/29 18:36:36 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestOverlapList.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestOverlapList.cc,v 1.3 2006/06/29 18:36:39 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestOvershootList.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestOvershootList.cc,v 1.3 2006/06/29 18:36:41 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestPoint.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestPoint.cc,v 1.3 2006/06/29 18:36:44 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestSegment.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestSegment.cc,v 1.11 2007/11/16 09:39:14 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestStreamLogger.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestStreamLogger.cc,v 1.3 2006/06/29 18:36:49 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestVolPoint.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestVolPoint.cc,v 1.3 2006/06/29 18:36:52 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeomTestVolume.cc

    r850 r921  
    2626//
    2727// $Id: G4GeomTestVolume.cc,v 1.6 2007/11/16 09:39:14 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4GeometryMessenger.cc

    r850 r921  
    2626//
    2727// $Id: G4GeometryMessenger.cc,v 1.5 2006/06/29 18:36:57 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/src/G4MultiNavigator.cc

    r831 r921  
    2525//
    2626//
    27 // $Id: G4MultiNavigator.cc,v 1.7 2007/11/02 13:48:43 japost Exp $
     27// $Id: G4MultiNavigator.cc,v 1.8 2008/10/24 14:00:03 gcosmo Exp $
    2828// GEANT4 tag $ Name:  $
    2929//
    … …  
    423423
    424424G4double G4MultiNavigator::ComputeSafety( const G4ThreeVector& position,
    425                                                 G4double       maxDistance)
     425                                          const G4double       maxDistance,
     426                                          const G4bool         state)
    426427{
    427428    // Recompute safety for the relevant point
    … …  
    434435    for( register int num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
    435436    {
    436        safety = (*pNavigatorIter)->ComputeSafety( position, maxDistance );
     437       safety = (*pNavigatorIter)->ComputeSafety( position, maxDistance, state);
    437438       if( safety < minSafety ) { minSafety = safety; }
    438439    }

  • trunk/source/geometry/navigation/src/G4Navigator.cc

    r831 r921  
    2525//
    2626//
    27 // $Id: G4Navigator.cc,v 1.37 2007/10/18 14:18:36 gcosmo Exp $
     27// $Id: G4Navigator.cc,v 1.38 2008/10/24 14:00:03 gcosmo Exp $
    2828// GEANT4 tag $ Name:  $
    2929//
    … …  
    12141214//
    12151215G4double G4Navigator::ComputeSafety( const G4ThreeVector &pGlobalpoint,
    1216                                      const G4double pMaxLength)
     1216                                     const G4double pMaxLength,
     1217                                     const G4bool keepState)
    12171218{
    12181219  G4double newSafety = 0.0;
    … …  
    12351236  }
    12361237#endif
     1238
     1239  if (keepState)  { SetSavedState(); }
    12371240
    12381241  G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2();
    … …  
    13221325  fPreviousSafety = newSafety;
    13231326
     1327  if (keepState)  { RestoreSavedState(); }
     1328
    13241329#ifdef G4DEBUG_NAVIGATION
    13251330  if( fVerbose > 1 )

  • trunk/source/geometry/navigation/src/G4NormalNavigation.cc

    r850 r921  
    2626//
    2727// $Id: G4NormalNavigation.cc,v 1.9 2007/05/11 13:43:59 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/src/G4ParameterisedNavigation.cc

    r850 r921  
    2626//
    2727// $Id: G4ParameterisedNavigation.cc,v 1.12 2007/11/09 16:06:02 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/src/G4PathFinder.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4PathFinder.cc,v 1.59 2008/04/29 15:32:54 gcosmo Exp $
     27// $Id: G4PathFinder.cc,v 1.61 2008/11/13 12:59:26 gcosmo Exp $
    2828// GEANT4 tag $ Name:  $
    2929//
    … …  
    581581  if( (!fNewTrack) && ( longMoveEnd && longMoveSaf ) )
    582582  { 
    583      G4ThreeVector  LastSafetyLocation;
    584        // Copy to keep last value - and restore
    585 
    586      LastSafetyLocation= fSafetyLocation;
    587 
    588583     // Recompute ComputeSafety for end position
    589584     //
    590585     revisedSafety= ComputeSafety(lastEndPosition);
    591 
    592      // Reset the state of last call to ComputeSafety
    593      //
    594      ComputeSafety( LastSafetyLocation );
    595586
    596587#ifdef G4DEBUG_PATHFINDER
    … …  
    757748 
    758749   std::vector<G4Navigator*>::iterator pNavigatorIter;
    759    pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator();
     750   pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
    760751
    761752   for( register G4int num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num )
    762753   {
    763       G4double safety = (*pNavigatorIter)->ComputeSafety( position );
     754      G4double safety = (*pNavigatorIter)->ComputeSafety( position,true );
    764755      if( safety < minSafety ) { minSafety = safety; }
    765756      fNewSafetyComputed[num]= safety;
    … …  
    11831174     for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
    11841175     {
    1185         safety= fpNavigator[numNav]->ComputeSafety( startPoint );
     1176        safety= fpNavigator[numNav]->ComputeSafety( startPoint, false );
    11861177        fPreSafetyValues[numNav]= safety;
    11871178        fCurrentPreStepSafety[numNav]= safety;

  • trunk/source/geometry/navigation/src/G4PropagatorInField.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4PropagatorInField.cc,v 1.43 2008/05/28 09:12:23 tnikitin Exp $
    28 // GEANT4 tag $Name: HEAD $
    2927//
    3028//
    … …  
    5048#include "G4VCurvedTrajectoryFilter.hh"
    5149#include "G4ChordFinder.hh"
     50#include "G4BrentLocator.hh"
    5251
    5352///////////////////////////////////////////////////////////////////////////
    … …  
    5655
    5756G4PropagatorInField::G4PropagatorInField( G4Navigator    *theNavigator,
    58                                           G4FieldManager *detectorFieldMgr )
     57                                          G4FieldManager *detectorFieldMgr,
     58                                          G4VIntersectionLocator *vLocator  )
    5959  : fDetectorFieldMgr(detectorFieldMgr),
    6060    fCurrentFieldMgr(detectorFieldMgr),
    … …  
    8383  fPreviousSafety= 0.0;
    8484  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    85   //
    86   fUseBrentLocator=true;
    87   // In case of too slow progress in finding Intersection Point
    88   // intermediates Points on the Track must be stored.
    89   // Initialise the array of Pointers [max_depth+1] to do this 
    90   G4ThreeVector zeroV(0.0,0.0,0.0);
    91   for (G4int idepth=0; idepth<max_depth+1; idepth++ )
    92   {
    93     ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
    94   }
    95   // Counter for Maximum Number Of Trial before Intersection Found
    96     maxNumberOfStepsForIntersection=0;
    97  // Counter for Number Of Calls to ReIntegrationEndPoint Method
    98     maxNumberOfCallsToReIntegration=0;
    99     maxNumberOfCallsToReIntegration_depth=0;
     85
     86  // Definding Intersection Locator and his parameters
     87  if(vLocator==0){
     88    fIntersectionLocator= new G4BrentLocator(theNavigator);
     89    fAllocatedLocator=true;
     90  }else{
     91    fIntersectionLocator=vLocator;
     92    fAllocatedLocator=false;
     93  }
     94  fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
     95  fIntersectionLocator->SetDeltaIntersectionFor(GetDeltaIntersection());
     96  fIntersectionLocator->SetChordFinderFor(GetChordFinder());
     97  fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
    10098}
    10199
    102100G4PropagatorInField::~G4PropagatorInField()
    103101{
    104   for ( G4int idepth=0; idepth<max_depth+1; idepth++)
    105   {
    106     delete ptrInterMedFT[idepth];
    107   }
    108   if(fVerboseLevel>0){
    109   G4cout<<"G4PropagatorInField::Location with Max Number of Steps="
    110        << maxNumberOfStepsForIntersection<<G4endl;
    111   G4cout<<"G4PropagatorInField::ReIntegrateEndPoint was called "<<
    112     maxNumberOfCallsToReIntegration<<" times and for depth algorithm "<<
    113     maxNumberOfCallsToReIntegration_depth<<" times"<<G4endl;
    114   }
     102  if(fAllocatedLocator)delete  fIntersectionLocator;
    115103}
    116104
    … …  
    162150  fSetFieldMgr= false;
    163151
    164   GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass); 
    165  
     152  GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass);
     153
     154 // Values for Intersection Locator has to be updated on each call
     155 // because the CurrentFieldManager changes
     156    fIntersectionLocator->SetChordFinderFor(GetChordFinder());
     157    fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
     158    fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
     159    fIntersectionLocator->SetDeltaIntersectionFor(GetDeltaIntersection());
     160
    166161  G4FieldTrack  CurrentState(pFieldTrack);
    167162  G4FieldTrack  OriginalState = CurrentState;
    … …  
    300295       //   of vol(A), if it exists. Start with point E as first "estimate".
    301296       G4bool recalculatedEndPt= false;
    302        G4bool found_intersection =
    303          LocateIntersectionPoint( SubStepStartState, CurrentState,
     297       
     298         G4bool found_intersection = fIntersectionLocator->
     299         EstimateIntersectionPoint( SubStepStartState, CurrentState,
    304300                                  InterSectionPointE, IntersectPointVelct_G,
    305                                   recalculatedEndPt);
    306        // G4cout<<"In Locate"<<recalculatedEndPt<<"  and V"<<IntersectPointVelct_G.GetPosition()<<G4endl;
     301                                  recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
    307302       intersects = intersects && found_intersection;
    308303       if( found_intersection ) {       
    … …  
    432427 
    433428  return TruePathLength;
    434 }
    435 
    436 // --------------------------------------------------------------------------
    437 // G4bool
    438 // G4PropagatorInField::LocateIntersectionPoint(
    439 //   const G4FieldTrack&       CurveStartPointVelocity,   //  A
    440 //   const G4FieldTrack&       CurveEndPointVelocity,     //  B
    441 //   const G4ThreeVector&      TrialPoint,                //  E
    442 //         G4FieldTrack&       IntersectedOrRecalculated  // Output
    443 //         G4bool&             recalculated)              // Out
    444 // --------------------------------------------------------------------------
    445 //
    446 // Function that returns the intersection of the true path with the surface
    447 // of the current volume (either the external one or the inner one with one
    448 // of the daughters
    449 //
    450 //     A = Initial point
    451 //     B = another point
    452 //
    453 // Both A and B are assumed to be on the true path.
    454 //
    455 //     E is the first point of intersection of the chord AB with
    456 //     a volume other than A  (on the surface of A or of a daughter)
    457 //
    458 // Convention of Use :
    459 //     i) If it returns "true", then IntersectionPointVelocity is set
    460 //       to the approximate intersection point.
    461 //    ii) If it returns "false", no intersection was found.
    462 //          The validity of IntersectedOrRecalculated depends on 'recalculated'
    463 //        a) if latter is false, then IntersectedOrRecalculated is invalid.
    464 //        b) if latter is true,  then IntersectedOrRecalculated is
    465 //             the new endpoint, due to a re-integration.
    466 // --------------------------------------------------------------------------
    467 
    468 G4bool
    469 G4PropagatorInField::LocateIntersectionPoint(
    470   const   G4FieldTrack&       CurveStartPointVelocity,   //  A
    471   const   G4FieldTrack&       CurveEndPointVelocity,     //  B
    472   const   G4ThreeVector&      TrialPoint,                //  E
    473           G4FieldTrack&       IntersectedOrRecalculatedFT, // Out: point found
    474           G4bool&             recalculatedEndPoint)        // Out:
    475 {
    476   // Find Intersection Point ( A, B, E )  of true path AB - start at E.
    477 
    478   G4bool found_approximate_intersection = false;
    479   G4bool there_is_no_intersection       = false;
    480  
    481   G4FieldTrack  CurrentA_PointVelocity = CurveStartPointVelocity;
    482   G4FieldTrack  CurrentB_PointVelocity = CurveEndPointVelocity;
    483   G4ThreeVector CurrentE_Point = TrialPoint;
    484   G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
    485   G4double    NewSafety= -0.0;
    486  
    487   G4bool final_section= true;  // Shows whether current section is last
    488                                // (i.e. B=full end)
    489   G4bool first_section=true;
    490   recalculatedEndPoint= false;
    491  
    492   G4bool restoredFullEndpoint= false;
    493 
    494   G4int substep_no = 0;
    495    
    496   // Limits for substep number
    497   //
    498   const G4int max_substeps=   10000;  // Test 120  (old value 100 )
    499   const G4int warn_substeps=   1000;  //      100 
    500 
    501   // Statistics for substeps
    502   //
    503   static G4int max_no_seen= -1;
    504   static G4int trigger_substepno_print= warn_substeps - 20 ;
    505  
    506   //-------------------------------------------------------------------------- 
    507   //  Algoritm for the case if progress in founding intersection is too slow.
    508   //  Process is defined too slow if after N=param_substeps advances on the
    509   //  path, it will be only 'fraction_done' of the total length.
    510   //  In this case the remaining length is divided in two half and
    511   //  the loop is restarted for each half. 
    512   //  If progress is still too slow, the division in two halfs continue
    513   //  until 'max_depth'.
    514   //--------------------------------------------------------------------------
    515   G4double count_did_len=0.;
    516   G4double count_all_len=0;
    517   G4int param_substeps=100;//Test value for the maximum number of substeps
    518   if(!fUseBrentLocator)  param_substeps=10;// Reduced value for the maximum number
    519 
    520   const G4double fraction_done=0.3;
    521 
    522   G4bool Second_half=false;      // First half or second half of divided step
    523 
    524   // We need to know this for the 'final_section':
    525   // real 'final_section' or first half 'final_section'
    526   // In algorithm it is considered that the 'Second_half' is true
    527   // and it becomes false only if we are in the first-half of level
    528   // depthness or if we are in the first section
    529 
    530   G4int depth=0; // Depth counts how many subdivisions of initial step made
    531 
    532 #ifdef G4DEBUG_FIELD
    533   static G4double tolerance= 1.0e-8;
    534   G4ThreeVector  StartPosition= CurveStartPointVelocity.GetPosition();
    535   if( (TrialPoint - StartPosition).mag() < tolerance * mm )
    536   {
    537      G4cerr << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
    538             << G4endl
    539             << "          Intermediate F point is on top of starting point A."
    540             << G4endl;
    541      G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
    542                  "IntersectionPointIsAtStart", JustWarning,
    543                  "Intersection point F is exactly at start point A." );
    544   }
    545 #endif
    546 
    547   // Intermediates Points on the Track = Subdivided Points must be stored.
    548   // Give the initial values to 'InterMedFt'
    549   // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
    550   //
    551   *ptrInterMedFT[0] = CurveEndPointVelocity;
    552   for (G4int idepth=1; idepth<max_depth+1; idepth++ )
    553   {
    554     *ptrInterMedFT[idepth]=CurveStartPointVelocity;
    555   }
    556 
    557   // 'SubStartPoint' is needed to calculate the length of the divided step
    558   //
    559   G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
    560    
    561   do
    562   {
    563     G4int substep_no_p = 0;
    564     G4bool sub_final_section = false; // the same as final_section,
    565                                       // but for 'sub_section'
    566     do // REPEAT param
    567     {
    568       G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition(); 
    569       G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
    570        
    571       // F = a point on true AB path close to point E
    572       // (the closest if possible)
    573       //
    574       if((!fUseBrentLocator)||(substep_no_p==0)){
    575        ApproxIntersecPointV = GetChordFinder()
    576                              ->ApproxCurvePointV( CurrentA_PointVelocity,
    577                                                   CurrentB_PointVelocity,
    578                                                   CurrentE_Point,
    579                                                   fEpsilonStep );
    580       //  The above method is the key & most intuitive part ...
    581       }
    582 #ifdef G4DEBUG_FIELD
    583       if( ApproxIntersecPointV.GetCurveLength() >
    584           CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
    585       {
    586         G4cerr << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
    587                << G4endl
    588                << "        Intermediate F point is more advanced than"
    589                << " endpoint B." << G4endl;
    590         G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
    591                     "IntermediatePointConfusion", FatalException,
    592                     "Intermediate F point is past end B point" );
    593       }
    594 #endif
    595 
    596       G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
    597       if(substep_no> maxNumberOfStepsForIntersection)maxNumberOfStepsForIntersection=substep_no; 
    598       // First check whether EF is small - then F is a good approx. point
    599       // Calculate the length and direction of the chord AF
    600       //
    601       G4ThreeVector  ChordEF_Vector = CurrentF_Point - CurrentE_Point;
    602 
    603       if ( ChordEF_Vector.mag2() <= sqr(GetDeltaIntersection()) )
    604       {
    605         found_approximate_intersection = true;
    606         // Create the "point" return value
    607         //
    608         IntersectedOrRecalculatedFT = ApproxIntersecPointV;
    609         IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
    610        
    611         // Note: in order to return a point on the boundary,
    612         //       we must return E. But it is F on the curve.
    613         //       So we must "cheat": we are using the position at point E
    614         //       and the velocity at point F !!!
    615         //
    616         // This must limit the length we can allow for displacement!
    617       }
    618       else  // E is NOT close enough to the curve (ie point F)
    619       {
    620         // Check whether any volumes are encountered by the chord AF
    621         // ---------------------------------------------------------
    622         // First relocate to restore any Voxel etc information
    623         // in the Navigator before calling ComputeStep()
    624         //
    625         fNavigator->LocateGlobalPointWithinVolume( Point_A );
    626 
    627         G4ThreeVector PointG;   // Candidate intersection point
    628         G4double stepLengthAF;
    629         G4bool Intersects_AF = IntersectChord( Point_A,   CurrentF_Point,
    630                                                NewSafety, stepLengthAF,
    631                                                PointG );
    632         if( Intersects_AF )
    633         {
    634           if(fUseBrentLocator){
    635            
    636             G4FieldTrack EndPoint=ApproxIntersecPointV;
    637             ApproxIntersecPointV= GetChordFinder()->ApproxCurvePointS(
    638                    CurrentA_PointVelocity,CurrentB_PointVelocity,
    639                    CurrentE_Point,CurrentF_Point,PointG,true,fEpsilonStep);
    640             CurrentB_PointVelocity =  EndPoint;
    641             CurrentE_Point = PointG;
    642           // By moving point B, must take care if current
    643           // AF has no intersection to try current FB!!
    644           //
    645           final_section= false;
    646 
    647           }
    648           else{
    649           // G is our new Candidate for the intersection point.
    650           // It replaces  "E" and we will repeat the test to see if
    651           // it is a good enough approximate point for us.
    652           //       B    <- F
    653           //       E    <- G
    654 
    655           CurrentB_PointVelocity = ApproxIntersecPointV;
    656           CurrentE_Point = PointG; 
    657      
    658           // By moving point B, must take care if current
    659           // AF has no intersection to try current FB!!
    660           //
    661           final_section= false;
    662           }
    663 #ifdef G4VERBOSE
    664           if( fVerboseLevel > 3 )
    665           {
    666             G4cout << "G4PiF::LI> Investigating intermediate point"
    667                    << " at s=" << ApproxIntersecPointV.GetCurveLength()
    668                    << " on way to full s="
    669                    << CurveEndPointVelocity.GetCurveLength() << G4endl;
    670           }
    671 #endif
    672         }
    673         else  // not Intersects_AF
    674         { 
    675           // In this case:
    676           // There is NO intersection of AF with a volume boundary.
    677           // We must continue the search in the segment FB!
    678           //
    679           fNavigator->LocateGlobalPointWithinVolume( CurrentF_Point );
    680 
    681           G4double stepLengthFB;
    682           G4ThreeVector PointH;
    683 
    684           // Check whether any volumes are encountered by the chord FB
    685           // ---------------------------------------------------------
    686 
    687           G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
    688                                                  NewSafety, stepLengthFB,
    689                                                  PointH );
    690           if( Intersects_FB )
    691           {
    692             if(fUseBrentLocator){
    693                CurrentA_PointVelocity = ApproxIntersecPointV;
    694                ApproxIntersecPointV= GetChordFinder()->ApproxCurvePointS(
    695                    CurrentA_PointVelocity,CurrentB_PointVelocity,
    696                    CurrentE_Point,Point_A,PointH,false,fEpsilonStep);
    697                CurrentE_Point = PointH;
    698            }
    699            else{
    700  
    701             // There is an intersection of FB with a volume boundary
    702             // H <- First Intersection of Chord FB
    703 
    704             // H is our new Candidate for the intersection point.
    705             // It replaces  "E" and we will repeat the test to see if
    706             // it is a good enough approximate point for us.
    707 
    708             // Note that F must be in volume volA  (the same as A)
    709             // (otherwise AF would meet a volume boundary!)
    710             //   A    <- F
    711             //   E    <- H
    712 
    713             CurrentA_PointVelocity = ApproxIntersecPointV;
    714             CurrentE_Point = PointH;
    715            }
    716           }
    717           else  // not Intersects_FB
    718           {
    719             // There is NO intersection of FB with a volume boundary
    720 
    721             if( final_section  )
    722             {
    723               // If B is the original endpoint, this means that whatever
    724               // volume(s) intersected the original chord, none touch the
    725               // smaller chords we have used.
    726               // The value of 'IntersectedOrRecalculatedFT' returned is
    727               // likely not valid
    728 
    729               // Check on real final_section or SubEndSection
    730               //
    731               if( ((Second_half)&&(depth==0)) || (first_section) )
    732               {
    733                 there_is_no_intersection = true;   // real final_section
    734               }
    735               else
    736               {
    737                 // end of subsection, not real final section
    738                 // exit from the and go to the depth-1 level
    739 
    740                 substep_no_p = param_substeps+2;  // exit from the loop
    741 
    742                 // but 'Second_half' is still true because we need to find
    743                 // the 'CurrentE_point' for the next loop
    744                 //
    745                 Second_half = true;
    746                 sub_final_section = true;
    747            
    748               }
    749             }
    750             else
    751             {
    752               // We must restore the original endpoint
    753 
    754               CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
    755               CurrentB_PointVelocity = CurveEndPointVelocity;
    756               restoredFullEndpoint = true;
    757             }
    758           } // Endif (Intersects_FB)
    759         } // Endif (Intersects_AF)
    760 
    761         // Ensure that the new endpoints are not further apart in space
    762         // than on the curve due to different errors in the integration
    763         //
    764         G4double linDistSq, curveDist;
    765         linDistSq = ( CurrentB_PointVelocity.GetPosition()
    766                     - CurrentA_PointVelocity.GetPosition() ).mag2();
    767         curveDist = CurrentB_PointVelocity.GetCurveLength()
    768                     - CurrentA_PointVelocity.GetCurveLength();
    769          if( curveDist*curveDist*(1+2*fEpsilonStep ) < linDistSq )
    770         {
    771           // Re-integrate to obtain a new B
    772           //
    773           G4FieldTrack newEndPointFT=
    774                   ReEstimateEndpoint( CurrentA_PointVelocity,
    775                                       CurrentB_PointVelocity,
    776                                       linDistSq,    // to avoid recalculation
    777                                       curveDist );
    778           G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
    779           CurrentB_PointVelocity = newEndPointFT;
    780           maxNumberOfCallsToReIntegration= maxNumberOfCallsToReIntegration+1;
    781           #ifdef G4DEBUG_FIELD
    782           G4cout<<"G4PIF::Call ReIntEnd1 linD="<<std::sqrt(linDistSq)<<" curve="<<curveDist<<" n="<<substep_no<<G4endl;
    783           G4cout<<"G4PIF::Call ReIntEnd2 IntersectAF="<< Intersects_AF<<" final_section="<<final_section<<G4endl;
    784           #endif
    785           if( (final_section)&&(Second_half)&&(depth==0) ) // real final section
    786           {
    787             recalculatedEndPoint = true;
    788             IntersectedOrRecalculatedFT = newEndPointFT;
    789               // So that we can return it, if it is the endpoint!
    790           }
    791          
    792         }
    793    
    794         if( curveDist < 0.0 )
    795         {
    796           G4cerr << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
    797                  << G4endl
    798                  << "        Error in advancing propagation." << G4endl;
    799           fVerboseLevel = 5; // Print out a maximum of information
    800           printStatus( CurrentA_PointVelocity,  CurrentB_PointVelocity,
    801                        -1.0, NewSafety,  substep_no, 0 );
    802           G4cerr << "        Point A (start) is " << CurrentA_PointVelocity
    803                  << G4endl;
    804           G4cerr << "        Point B (end)   is " << CurrentB_PointVelocity
    805                  << G4endl;
    806           G4cerr << "        Curve distance is " << curveDist << G4endl;
    807           G4cerr << G4endl
    808                  << "The final curve point is not further along"
    809                  << " than the original!" << G4endl;
    810           if( recalculatedEndPoint )
    811           {
    812             G4cerr << "Recalculation of EndPoint was called with fEpsStep= "
    813                    << fEpsilonStep << G4endl;
    814           }
    815           G4cerr.precision(20);
    816           G4cerr << " Point A (Curve start)   is " << CurveStartPointVelocity
    817                  << G4endl;
    818           G4cerr << " Point B (Curve   end)   is " << CurveEndPointVelocity
    819                  << G4endl;
    820           G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity
    821                  << G4endl;
    822           G4cerr << " Point B (Current end)   is " << CurrentB_PointVelocity
    823                  << G4endl;
    824           G4cerr << " Point S (Sub start)     is " << SubStart_PointVelocity
    825                  << G4endl;
    826           G4cerr << " Point E (Trial Point)   is " << CurrentE_Point
    827                  << G4endl;
    828           G4cerr << " Point F (Intersection)  is " << ApproxIntersecPointV
    829                  << G4endl;
    830           G4cerr << "        LocateIntersection parameters are : Substep no= "
    831                  << substep_no << G4endl;
    832           G4cerr << "        Substep depth no= "<< substep_no_p  << " Depth= "
    833                  << depth << G4endl;
    834           G4cerr << "        did_len= "<< count_did_len  << " all_len= "
    835                  << count_all_len << G4endl;
    836           G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
    837                       "FatalError", FatalException,
    838                       "Error in advancing propagation.");
    839         }
    840        
    841         if(restoredFullEndpoint)
    842         {
    843           final_section = restoredFullEndpoint;
    844           restoredFullEndpoint = false;
    845         }
    846       } // EndIf ( E is close enough to the curve, ie point F. )
    847         // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
    848 
    849 #ifdef G4DEBUG_LOCATE_INTERSECTION 
    850       if( substep_no >= trigger_substepno_print )
    851       {
    852         G4cout << "Difficulty in converging in "
    853                << "G4PropagatorInField::LocateIntersectionPoint():"
    854                << G4endl
    855                << "    Substep no = " << substep_no << G4endl;
    856         if( substep_no == trigger_substepno_print )
    857         {
    858           printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
    859                        -1.0, NewSafety, 0, 0);
    860         }
    861         G4cout << " State of point A: ";
    862         printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
    863                      -1.0, NewSafety, substep_no-1, 0);
    864         G4cout << " State of point B: ";
    865         printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
    866                      -1.0, NewSafety, substep_no, 0);
    867       }
    868 #endif
    869 
    870       substep_no++;
    871       substep_no_p++;
    872 
    873     } while (  ( ! found_approximate_intersection )
    874             && ( ! there_is_no_intersection )     
    875             && ( substep_no_p <= param_substeps) );  // UNTIL found or
    876                                                      // failed param substep
    877     first_section = false;
    878 
    879     if( (!found_approximate_intersection) && (!there_is_no_intersection) )
    880     {
    881       G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
    882                        - SubStart_PointVelocity.GetCurveLength());
    883       G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
    884                        - SubStart_PointVelocity.GetCurveLength());
    885       count_did_len=did_len;
    886       count_all_len=all_len;   
    887       G4double stepLengthAB;
    888       G4ThreeVector PointGe;
    889 
    890       // Check if progress is too slow and if it possible to go deeper,
    891       // then halve the step if so
    892       //
    893       if( ( ( did_len )<fraction_done*all_len)
    894          && (depth<max_depth) && (!sub_final_section) )
    895       {
    896 
    897         Second_half=false;
    898         depth++;
    899 
    900         G4double Sub_len = (all_len-did_len)/(2.);
    901         G4FieldTrack start = CurrentA_PointVelocity;
    902         G4MagInt_Driver* integrDriver=GetChordFinder()->GetIntegrationDriver();
    903         integrDriver->AccurateAdvance(start, Sub_len, fEpsilonStep);
    904         *ptrInterMedFT[depth] = start;
    905         CurrentB_PointVelocity = *ptrInterMedFT[depth];
    906  
    907         // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
    908         //
    909         SubStart_PointVelocity = CurrentA_PointVelocity;
    910 
    911         // Find new trial intersection point needed at start of the loop
    912         //
    913         G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
    914         G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();   
    915      
    916         fNavigator->LocateGlobalPointWithinVolume(Point_A);
    917         G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
    918                                               NewSafety, stepLengthAB, PointGe);
    919         if(Intersects_AB)
    920         {
    921           CurrentE_Point = PointGe;
    922         }
    923         else
    924         {
    925           // No intersection found for first part of curve
    926           // (CurrentA,InterMedPoint[depth]). Go to the second part
    927           //
    928           Second_half = true;
    929         }
    930       } // if did_len
    931 
    932       if( (Second_half)&&(depth!=0) )
    933       {
    934         // Second part of curve (InterMed[depth],Intermed[depth-1])                       )
    935         // On the depth-1 level normally we are on the 'second_half'
    936 
    937         Second_half = true;
    938 
    939         //  Find new trial intersection point needed at start of the loop
    940         //
    941         SubStart_PointVelocity = *ptrInterMedFT[depth];
    942         CurrentA_PointVelocity = *ptrInterMedFT[depth];
    943         CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
    944         // Ensure that the new endpoints are not further apart in space
    945         // than on the curve due to different errors in the integration
    946         //
    947         G4double linDistSq, curveDist;
    948         linDistSq = ( CurrentB_PointVelocity.GetPosition()
    949                     - CurrentA_PointVelocity.GetPosition() ).mag2();
    950         curveDist = CurrentB_PointVelocity.GetCurveLength()
    951                     - CurrentA_PointVelocity.GetCurveLength();
    952         if( curveDist*curveDist*(1+2*fEpsilonStep ) < linDistSq )
    953         {
    954           // Re-integrate to obtain a new B
    955           //
    956           G4FieldTrack newEndPointFT=
    957                   ReEstimateEndpoint( CurrentA_PointVelocity,
    958                                       CurrentB_PointVelocity,
    959                                       linDistSq,    // to avoid recalculation
    960                                       curveDist );
    961           G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
    962           CurrentB_PointVelocity = newEndPointFT;
    963           maxNumberOfCallsToReIntegration_depth= maxNumberOfCallsToReIntegration_depth+1;
    964         }
    965 
    966         G4ThreeVector Point_A    = CurrentA_PointVelocity.GetPosition();
    967         G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();   
    968         fNavigator->LocateGlobalPointWithinVolume(Point_A);
    969         G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
    970                                               stepLengthAB, PointGe);
    971         if(Intersects_AB)
    972         {
    973           CurrentE_Point = PointGe;
    974         }
    975         else
    976         {
    977           final_section = true;
    978         }
    979         depth--;
    980       }
    981     }  // if(!found_aproximate_intersection)
    982 
    983   } while ( ( ! found_approximate_intersection )
    984             && ( ! there_is_no_intersection )     
    985             && ( substep_no <= max_substeps) ); // UNTIL found or failed
    986 
    987   if( substep_no > max_no_seen )
    988   {
    989     max_no_seen = substep_no;
    990     if( max_no_seen > warn_substeps )
    991     {
    992       trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
    993     }
    994   }
    995 
    996   if(  ( substep_no >= max_substeps)
    997       && !there_is_no_intersection
    998       && !found_approximate_intersection )
    999   {
    1000     G4cerr << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
    1001            << G4endl
    1002            << "          Convergence is requiring too many substeps: "
    1003            << substep_no << G4endl;
    1004     G4cerr << "          Abandoning effort to intersect. " << G4endl;
    1005     G4cerr << "          Information on start & current step follows in cout."
    1006            << G4endl;
    1007     G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
    1008            << G4endl
    1009            << "          Convergence is requiring too many substeps: "
    1010            << substep_no << G4endl;
    1011     G4cout << "          Found intersection = "
    1012            << found_approximate_intersection << G4endl
    1013            << "          Intersection exists = "
    1014            << !there_is_no_intersection << G4endl;
    1015     G4cout << "          Start and Endpoint of Requested Step:" << G4endl;
    1016     printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
    1017                  -1.0, NewSafety, 0, 0);
    1018     G4cout << G4endl;
    1019     G4cout << "          'Bracketing' starting and endpoint of current Sub-Step"
    1020            << G4endl;
    1021     printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
    1022                  -1.0, NewSafety, substep_no-1, 0);
    1023     printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
    1024                  -1.0, NewSafety, substep_no, 0);
    1025     G4cout << G4endl;
    1026  
    1027 #ifdef FUTURE_CORRECTION
    1028     // Attempt to correct the results of the method // FIX - TODO
    1029 
    1030     if ( ! found_approximate_intersection )
    1031     {
    1032       recalculatedEndPoint = true;
    1033       // Return the further valid intersection point -- potentially A ??
    1034       // JA/19 Jan 2006
    1035       IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
    1036 
    1037       G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
    1038              << G4endl
    1039              << "          Did not convergence after " << substep_no
    1040              << " substeps." << G4endl;
    1041       G4cout << "          The endpoint was adjused to pointA resulting"
    1042              << G4endl
    1043              << "          from the last substep: " << CurrentA_PointVelocity
    1044              << G4endl;
    1045     }
    1046 #endif
    1047 
    1048     G4cout.precision( 10 );
    1049     G4double done_len = CurrentA_PointVelocity.GetCurveLength();
    1050     G4double full_len = CurveEndPointVelocity.GetCurveLength();
    1051     G4cout << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
    1052            << G4endl
    1053            << "        Undertaken only length: " << done_len
    1054            << " out of " << full_len << " required." << G4endl;
    1055     G4cout << "        Remaining length = " << full_len - done_len << G4endl;
    1056 
    1057     G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
    1058                 "UnableToLocateIntersection", FatalException,
    1059                 "Too many substeps while trying to locate intersection.");
    1060   }
    1061   else if( substep_no >= warn_substeps )
    1062   { 
    1063     int oldprc= G4cout.precision( 10 );
    1064     G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
    1065            << G4endl
    1066            << "          Undertaken length: " 
    1067            << CurrentB_PointVelocity.GetCurveLength();
    1068     G4cout << " - Needed: "  << substep_no << " substeps." << G4endl
    1069            << "          Warning level = " << warn_substeps
    1070            << " and maximum substeps = " << max_substeps << G4endl;
    1071     G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
    1072                 "DifficultyToLocateIntersection", JustWarning,
    1073                 "Many substeps while trying to locate intersection.");
    1074     G4cout.precision( oldprc );
    1075   }
    1076  
    1077   return  !there_is_no_intersection; //  Success or failure
    1078429}
    1079430
    … …  
    1210561}
    1211562
    1212 G4bool
    1213 G4PropagatorInField::IntersectChord( G4ThreeVector  StartPointA,
    1214                                      G4ThreeVector  EndPointB,
    1215                                      G4double      &NewSafety,
    1216                                      G4double      &LinearStepLength,
    1217                                      G4ThreeVector &IntersectionPoint
    1218                                    )
    1219 {
    1220     // Calculate the direction and length of the chord AB
    1221     G4ThreeVector  ChordAB_Vector = EndPointB - StartPointA;
    1222     G4double       ChordAB_Length = ChordAB_Vector.mag();  // Magnitude (norm)
    1223     G4ThreeVector  ChordAB_Dir =    ChordAB_Vector.unit();
    1224     G4bool intersects;
    1225 
    1226     G4ThreeVector OriginShift = StartPointA - fPreviousSftOrigin ;
    1227     G4double      MagSqShift  = OriginShift.mag2() ;
    1228     G4double      currentSafety;
    1229     G4bool        doCallNav= false;
    1230 
    1231     if( MagSqShift >= sqr(fPreviousSafety) )
    1232     {
    1233         currentSafety = 0.0 ;
    1234     }else{
    1235         currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
    1236     }
    1237 
    1238     if( fUseSafetyForOptimisation && (ChordAB_Length <= currentSafety) )
    1239     {
    1240        // The Step is guaranteed to be taken
    1241 
    1242        LinearStepLength = ChordAB_Length;
    1243        intersects = false;
    1244 
    1245        NewSafety= currentSafety;
    1246 
    1247 #if 0
    1248        G4cout << " G4PropagatorInField does not call Navigator::ComputeStep " << G4endl ;
    1249        G4cout << "    step= " << LinearStepLength << " safety= " << NewSafety << G4endl;
    1250        G4cout << "    safety: Origin = " << fPreviousSftOrigin << " val= " << fPreviousSafety << G4endl;
    1251 #endif
    1252     }
    1253     else
    1254     {
    1255        doCallNav= true;
    1256        // Check whether any volumes are encountered by the chord AB
    1257 
    1258        // G4cout << " G4PropagatorInField calling Navigator::ComputeStep " << G4endl ;
    1259 
    1260        LinearStepLength =
    1261         fNavigator->ComputeStep( StartPointA, ChordAB_Dir,
    1262                                  ChordAB_Length, NewSafety );
    1263        intersects = (LinearStepLength <= ChordAB_Length);
    1264        // G4Navigator contracts to return k_infinity if len==asked
    1265        // and it did not find a surface boundary at that length
    1266        LinearStepLength = std::min( LinearStepLength, ChordAB_Length);
    1267 
    1268        // G4cout << " G4PiF got step= " << LinearStepLength << " safety= " << NewSafety << G4endl;
    1269 
    1270        // Save the last calculated safety!
    1271        fPreviousSftOrigin = StartPointA;
    1272        fPreviousSafety= NewSafety;
    1273 
    1274        if( intersects ){
    1275           // Intersection Point of chord AB and either volume A's surface
    1276           //                                or a daughter volume's surface ..
    1277           IntersectionPoint = StartPointA + LinearStepLength * ChordAB_Dir;
    1278        }
    1279     }
    1280 
    1281 #ifdef DEBUG_INTERSECTS_CHORD
    1282     // printIntersection(
    1283     // StartPointA, EndPointB, LinearStepLength, IntersectionPoint, NewSafety
    1284 
    1285     G4cout << " G4PropagatorInField::IntersectChord reports " << G4endl;
    1286     G4cout << " PiF-IC> "
    1287            << "Start="  << std::setw(12) << StartPointA       << " "
    1288            << "End= "   << std::setw(8) << EndPointB         << " "
    1289            << "StepIn=" << std::setw(8) << LinearStepLength  << " "
    1290            << "NewSft=" << std::setw(8) << NewSafety << " "
    1291            << "CallNav=" << doCallNav      << "  "
    1292            << "Intersects " << intersects     << "  ";
    1293     if( intersects )
    1294       G4cout << "IntrPt=" << std::setw(8) << IntersectionPoint << " " ;
    1295     G4cout << G4endl;
    1296 #endif
    1297 
    1298     return intersects;
    1299 }
    1300 
    1301 // --------------------- oooo000000000000oooo ----------------------------
    1302 
    1303 G4FieldTrack G4PropagatorInField::
    1304 ReEstimateEndpoint( const G4FieldTrack &CurrentStateA, 
    1305                     const G4FieldTrack &EstimatedEndStateB,
    1306                           G4double      linearDistSq,
    1307                           G4double      curveDist
    1308                   )
    1309 {
    1310   // G4double checkCurveDist= EstimatedEndStateB.GetCurveLength()
    1311   //   - CurrentStateA.GetCurveLength();
    1312   // G4double checkLinDistSq= (EstimatedEndStateB.GetPosition()
    1313   //                    - CurrentStateA.GetPosition() ).mag2();
    1314 
    1315   G4FieldTrack newEndPoint( CurrentStateA );
    1316   G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
    1317 
    1318   G4FieldTrack retEndPoint( CurrentStateA );
    1319   G4bool goodAdvance;
    1320   G4int  itrial=0;
    1321   const G4int no_trials= 20;
    1322 
    1323   G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
    1324   do
    1325   {
    1326      G4double currentCurveLen= newEndPoint.GetCurveLength();
    1327      G4double advanceLength= endCurveLen - currentCurveLen ;
    1328      if (std::abs(advanceLength)<kCarTolerance)
    1329      {
    1330        advanceLength=(EstimatedEndStateB.GetPosition()
    1331                      -newEndPoint.GetPosition()).mag();
    1332      }
    1333      goodAdvance=
    1334        integrDriver->AccurateAdvance(newEndPoint, advanceLength, fEpsilonStep);
    1335      //              ***************
    1336   }
    1337   while( !goodAdvance && (++itrial < no_trials) );
    1338 
    1339   if( goodAdvance )
    1340   {
    1341     retEndPoint= newEndPoint;
    1342   }
    1343   else
    1344   {
    1345     retEndPoint= EstimatedEndStateB; // Could not improve without major work !!
    1346   }
    1347 
    1348   //  All the work is done
    1349   //   below are some diagnostics only -- before the return!
    1350   //
    1351   static const G4String MethodName("G4PropagatorInField::ReEstimateEndpoint");
    1352 
    1353 #ifdef G4VERBOSE
    1354   G4int  latest_good_trials=0;
    1355   if( itrial > 1)
    1356   {
    1357     if( fVerboseLevel > 0 )
    1358     {
    1359       G4cout << MethodName << " called - goodAdv= " << goodAdvance
    1360              << " trials = " << itrial
    1361              << " previous good= " << latest_good_trials
    1362              << G4endl;
    1363     }
    1364     latest_good_trials=0;
    1365   }
    1366   else
    1367   {   
    1368     latest_good_trials++;
    1369   }
    1370 #endif
    1371 
    1372 #ifdef G4DEBUG_FIELD
    1373   G4double lengthDone = newEndPoint.GetCurveLength()
    1374                       - CurrentStateA.GetCurveLength();
    1375   if( !goodAdvance )
    1376   {
    1377     if( fVerboseLevel >= 3 )
    1378     {
    1379       G4cout << MethodName << "> AccurateAdvance failed " ;
    1380       G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
    1381       G4cout << " It went only " << lengthDone << " instead of " << curveDist
    1382              << " -- a difference of " << curveDist - lengthDone  << G4endl;
    1383       G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
    1384              << G4endl;
    1385     }
    1386   }
    1387 
    1388   static G4int noInaccuracyWarnings = 0;
    1389   G4int maxNoWarnings = 10;
    1390   if (  (noInaccuracyWarnings < maxNoWarnings )
    1391        || (fVerboseLevel > 1) )
    1392     {
    1393       G4cerr << "G4PropagatorInField::LocateIntersectionPoint():"
    1394              << G4endl
    1395              << " Warning: Integration inaccuracy requires"
    1396              <<   " an adjustment in the step's endpoint."  << G4endl
    1397              << "   Two mid-points are further apart than their"
    1398              <<   " curve length difference"                << G4endl
    1399              << "   Dist = "       << std::sqrt(linearDistSq)
    1400              << " curve length = " << curveDist             << G4endl;
    1401       G4cerr << " Correction applied is "
    1402              << (newEndPoint.GetPosition()-EstimatedEndStateB.GetPosition()).mag()
    1403              << G4endl;
    1404     }
    1405 #else
    1406   // Statistics on the RMS value of the corrections
    1407 
    1408   static G4int    noCorrections=0;
    1409   static G4double sumCorrectionsSq = 0;
    1410   noCorrections++;
    1411   if( goodAdvance )
    1412   {
    1413     sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
    1414                          newEndPoint.GetPosition()).mag2();
    1415   }
    1416   linearDistSq -= curveDist; // To use linearDistSq ... !
    1417 #endif
    1418 
    1419   return retEndPoint;
    1420 }
    1421 
    1422563// Access the points which have passed through the filter. The
    1423564// points are stored as ThreeVectors for the initial impelmentation

  • trunk/source/geometry/navigation/src/G4ReplicaNavigation.cc

    r850 r921  
    2626//
    2727// $Id: G4ReplicaNavigation.cc,v 1.19 2008/04/28 15:39:55 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/src/G4SafetyHelper.cc

    r831 r921  
    2424// ********************************************************************
    2525//
    26 // $Id: G4SafetyHelper.cc,v 1.15 2007/11/14 10:04:21 gcosmo Exp $
     26// $Id: G4SafetyHelper.cc,v 1.16 2008/10/24 14:00:03 gcosmo Exp $
    2727// GEANT4 tag $ Name:  $
    2828//
    … …  
    128128    {
    129129      // Safety for mass geometry
    130       fLastSafety = fpMassNavigator->ComputeSafety(position);
     130      fLastSafety = fpMassNavigator->ComputeSafety(position,true);
    131131    }
    132132    else
    133133    {
    134134      // Safety for all geometries
    135       fLastSafety = fpPathFinder->ComputeSafety( position );
     135      fLastSafety = fpPathFinder->ComputeSafety(position);
    136136    }
    137137    newSafety = fLastSafety;

  • trunk/source/geometry/navigation/src/G4TransportationManager.cc

    r850 r921  
    2626//
    2727// $Id: G4TransportationManager.cc,v 1.15 2007/04/12 11:51:48 vnivanch Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/src/G4VoxelNavigation.cc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4VoxelNavigation.cc,v 1.7 2007/05/11 13:43:59 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4VoxelNavigation.cc,v 1.9 2008/11/14 18:26:35 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    254254                G4String solidResponse = "-kInside-";
    255255                if (insideIntPt == kOutside)
    256                   solidResponse = "-kOutside-";
     256                  { solidResponse = "-kOutside-"; }
    257257                else if (insideIntPt == kSurface)
    258                   solidResponse = "-kSurface-";
     258                  { solidResponse = "-kSurface-"; }
    259259                if( fVerbose == 1 )
    260260                {
    … …  
    265265                         << "    For point p: " << intersectionPoint
    266266                         << ", considered as 'intersection' point." << G4endl;
     267                }
     268                G4double safetyIn= -1, safetyOut= -1;  //  Set to invalid values
     269                G4double newDistIn= -1,  newDistOut= -1;
     270                if( insideIntPt != kInside )
     271                {
     272                  safetyIn= sampleSolid->DistanceToIn(intersectionPoint);
     273                  newDistIn= sampleSolid->DistanceToIn(intersectionPoint,
     274                                                       sampleDirection);
     275                }
     276                if( insideIntPt != kOutside )
     277                {
     278                  safetyOut= sampleSolid->DistanceToOut(intersectionPoint);
     279                  newDistOut= sampleSolid->DistanceToOut(intersectionPoint,
     280                                                         sampleDirection);
    267281                }
    268282                if( insideIntPt != kSurface )
    … …  
    277291                         << " for this point !" << G4endl;
    278292                  G4cout << "          Point = " << intersectionPoint << G4endl;
     293                  G4cout << "          Safety values: " << G4endl;
    279294                  if ( insideIntPt != kInside )
    280                     G4cout << "        DistanceToIn(p) = "
    281                            << sampleSolid->DistanceToIn(intersectionPoint)
     295                  {
     296                    G4cout << "          DistanceToIn(p)  = " << safetyIn
    282297                           << G4endl;
    283                   if ( insideIntPt != kOutside )
    284                     G4cout << "        DistanceToOut(p) = "
    285                            << sampleSolid->DistanceToOut(intersectionPoint)
     298                  }
     299                  if ( insideIntPt != kOutside )
     300                  {
     301                    G4cout << "          DistanceToOut(p) = " << safetyOut
    286302                           << G4endl;
     303                  }
    287304                  G4Exception("G4VoxelNavigation::ComputeStep()",
    288305                              "InaccurateDistanceToIn", JustWarning,
    289                               "Navigator gets conflicting response from Solid.");
     306                              "Conflicting response from Solid.");
    290307                  G4cout.precision(oldcoutPrec);
     308                }
     309                else
     310                { 
     311                  // If it is on the surface, *ensure* that either DistanceToIn
     312                  // or DistanceToOut returns a finite value ( >= Tolerance).
     313                  //
     314                  if( std::max( newDistIn, newDistOut ) <= kCarTolerance )
     315                  {
     316                    G4cout << "ERROR - G4VoxelNavigation::ComputeStep()"
     317                       << G4endl
     318                       << "  Identified point for which the solid "
     319                       << sampleSolid->GetName() << G4endl
     320                       << "  has MAJOR problem:  " << G4endl
     321                       << "  --> Both DistanceToIn(p,v) and DistanceToOut(p,v) "
     322                       << "return Zero, an equivalent value or negative value."
     323                       << G4endl;
     324                    G4cout << "    Solid: " << sampleSolid << G4endl;
     325                    G4cout << "    Point p= " << intersectionPoint << G4endl;
     326                    G4cout << "    Direction v= " << sampleDirection << G4endl;
     327                    G4cout << "    DistanceToIn(p,v)     = " << newDistIn
     328                           << G4endl;
     329                    G4cout << "    DistanceToOut(p,v,..) = " << newDistOut
     330                           << G4endl;
     331                    G4cout << "    Safety values: " << G4endl;
     332                    G4cout << "      DistanceToIn(p)  = " << safetyIn
     333                           << G4endl;
     334                    G4cout << "      DistanceToOut(p) = " << safetyOut
     335                           << G4endl;
     336                    G4Exception("G4VoxelNavigation::ComputeStep()",
     337                              "DistanceToInAndOutAreZero", FatalException,
     338                              "Zero from both Solid DistanceIn and Out(p,v).");
     339                  }
    291340                }
    292341              }
Note: See TracChangeset for help on using the changeset viewer.