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
Files:
57 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/navigation/History

    r850 r921  
    1 $Id: History,v 1.119 2008/04/29 15:33:05 gcosmo Exp $
     1$Id: History,v 1.125 2008/11/14 18:26:53 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20November 14th, 2008 - T.Nikitina, J.Apostolakis (geomnav-V09-01-09)
     21-----------------------------------------------
     22- Introduced first implementation of new optional method in locator classes
     23  AdjustementOfFoundIntersection() using surface-normal of the intersecting
     24  solid to boost accuracy. Added the optional call to the new method in each
     25  concrete locator.
     26- Removed unnecessary accessors for Brent locator in G4PropagatorInField.
     27- G4VoxelNavigation: implemented additional check when running in "check"
     28  mode; if it is on the surface, ensure that it can move on next step;
     29  either DistanceToIn(p,v) or DistanceToOut(p,v) should return a finite
     30  value greater than the tolerance.
     31
     32November 10th, 2008 - G.Cosmo (geomnav-V09-01-08)
     33-----------------------------
     34- G4PathFinder: cleared unecessary calls to ComputeSafety() in ReLocate().
     35
     36October 28th, 2008 - T.Nikitina (geomnav-V09-01-07)
     37-------------------------------
     38- Moved method LocateIntersectionPoint() in G4PropagatorInField to a separate
     39  class G4VIntersectionLocator, now allowing to use different location
     40  algorithms: Brent, MultiLevel, Simple.
     41- New classes: G4VIntersectionLocator, G4SimpleLocator, G4BrentLocator and
     42  G4MultiLevelLocator. 
     43- Coworks with tag "field-V09-01-03".
     44
     45October 10th, 2008 - G.Cosmo (geomnav-V09-01-06)
     46----------------------------
     47- Introduced optional Boolean argument in G4Navigator::ComputeSafety() to
     48  allow for computation of safety without modifying the state restoring of
     49  the navigator.
     50- Modified accordingly the following classes, for calls to ComputeSafety():
     51  G4SafetyHelper, G4PathFinder (now calling ComputeSafety() with TRUE
     52  argument to preserve navigator's state), G4MultiNavigator and
     53  G4ErrorPropagationNavigator.
     54
     55May 5th, 2008 - T.Nikitina (geomnav-V09-01-05)
     56--------------------------
     57- Added Brent method for LocateIntersectionPoint() in G4PropagatorInField.
     58  The Brent method is now used as default and can be switched off through
     59  call to the proper accessor function SetBrentMethod().
     60- Requires related update to G4ChordFinder in geometry/magneticfield module
     61  included in tag "field-V09-01-02".
    1962
    2063April 29th, 2008 - M.Asai (geomnav-V09-01-04)

  • trunk/source/geometry/navigation/include/G4AuxiliaryNavServices.hh

    r850 r921  
    2626//
    2727// $Id: G4AuxiliaryNavServices.hh,v 1.4 2007/05/22 07:48:08 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4AuxiliaryNavServices.icc

    r850 r921  
    2626//
    2727// $Id: G4AuxiliaryNavServices.icc,v 1.4 2007/05/22 07:48:08 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4DrawVoxels.hh

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

  • trunk/source/geometry/navigation/include/G4ErrorPropagationNavigator.hh

    r850 r921  
    2525//
    2626//
    27 // $Id: G4ErrorPropagationNavigator.hh,v 1.1 2007/05/16 12:49:18 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4ErrorPropagationNavigator.hh,v 1.2 2008/10/24 14:00:03 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    6464
    6565    G4double ComputeSafety(const G4ThreeVector &globalpoint,
    66                            const G4double pProposedMaxLength = DBL_MAX);
     66                           const G4double pProposedMaxLength = DBL_MAX,
     67                           const G4bool keepState = false);
    6768      // Calls the navigation in the detector geometry and then checks
    6869      // if the distance to surface is smaller than the proposed safety

  • trunk/source/geometry/navigation/include/G4GeomTestErrorList.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestLogger.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestOverlapList.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestOvershootList.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestPoint.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestSegment.hh

    r850 r921  
    2626//
    2727// $Id: G4GeomTestSegment.hh,v 1.4 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/include/G4GeomTestStreamLogger.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestVolPoint.hh

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

  • trunk/source/geometry/navigation/include/G4GeomTestVolume.hh

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

  • trunk/source/geometry/navigation/include/G4GeometryMessenger.hh

    r850 r921  
    2626//
    2727// $Id: G4GeometryMessenger.hh,v 1.4 2006/06/29 18:35:59 gunter Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// --------------------------------------------------------------------

  • trunk/source/geometry/navigation/include/G4MultiNavigator.hh

    r850 r921  
    2525//
    2626//
    27 // $Id: G4MultiNavigator.hh,v 1.4 2007/05/21 15:36:25 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4MultiNavigator.hh,v 1.5 2008/10/24 14:00:03 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    113113
    114114  G4double ComputeSafety(const G4ThreeVector &globalpoint,
    115                          const G4double pProposedMaxLength = DBL_MAX);
     115                         const G4double pProposedMaxLength = DBL_MAX,
     116                         const G4bool keepState = false);
    116117    // Calculate the isotropic distance to the nearest boundary
    117118    // in any geometry from the specified point in the global coordinate

  • trunk/source/geometry/navigation/include/G4Navigator.hh

    r850 r921  
    2525//
    2626//
    27 // $Id: G4Navigator.hh,v 1.26 2007/10/18 14:18:36 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4Navigator.hh,v 1.27 2008/10/24 14:00:03 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    181181
    182182  virtual G4double ComputeSafety(const G4ThreeVector &globalpoint,
    183                                  const G4double pProposedMaxLength = DBL_MAX);
     183                                 const G4double pProposedMaxLength = DBL_MAX,
     184                                 const G4bool keepState = false);
    184185    // Calculate the isotropic distance to the nearest boundary from the
    185186    // specified point in the global coordinate system.

  • trunk/source/geometry/navigation/include/G4Navigator.icc

    r850 r921  
    2626//
    2727// $Id: G4Navigator.icc,v 1.15 2007/10/18 14:18:36 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4NormalNavigation.hh

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

  • trunk/source/geometry/navigation/include/G4NormalNavigation.icc

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

  • trunk/source/geometry/navigation/include/G4ParameterisedNavigation.hh

    r850 r921  
    2626//
    2727// $Id: G4ParameterisedNavigation.hh,v 1.6 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/include/G4ParameterisedNavigation.icc

    r850 r921  
    2626//
    2727// $Id: G4ParameterisedNavigation.icc,v 1.7 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/include/G4PathFinder.hh

    r850 r921  
    2525//
    2626// $Id: G4PathFinder.hh,v 1.34 2007/11/02 12:28:31 japost Exp $
    27 // GEANT4 tag $Name: HEAD $
     27// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2828//
    2929// class G4PathFinder

  • trunk/source/geometry/navigation/include/G4PhantomParameterisation.hh

    r850 r921  
    2626//
    2727// $Id: G4PhantomParameterisation.hh,v 1.4 2008/01/22 15:02:36 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4PhantomParameterisation.icc

    r850 r921  
    2626//
    2727// $Id: G4PhantomParameterisation.icc,v 1.1 2007/10/17 19:13:58 arce Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//--------------------------------------------------------------------

  • trunk/source/geometry/navigation/include/G4PropagatorInField.hh

    r850 r921  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PropagatorInField.hh,v 1.14 2008/05/28 09:11:59 tnikitin Exp $
    28 // GEANT4 tag $Name: HEAD $
    2926//
    30 // class G4PropagatorInField
    31 //
    32 // Class description:
     27// $Id: G4PropagatorInField.hh,v 1.17 2008/11/13 14:28:56 tnikitin Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
     29//
     30//
     31// Class G4PropagatorInField
     32//
     33// class description:
    3334//
    3435// This class performs the navigation/propagation of a particle/track
    … …  
    3637// For the calculation of the path, it relies on the class G4ChordFinder.
    3738//
    38 // Key Method:
    39 //              ComputeStep(..)
     39// Key Method: ComputeStep(..)
     40
    4041// History:
    4142// -------
    … …  
    5455#include "G4FieldTrack.hh"
    5556#include "G4FieldManager.hh"
     57#include "G4VIntersectionLocator.hh"
     58
    5659class G4ChordFinder;
    5760
    … …  
    6669
    6770   G4PropagatorInField( G4Navigator    *theNavigator,
    68                         G4FieldManager *detectorFieldMgr );
     71                        G4FieldManager *detectorFieldMgr,
     72                        G4VIntersectionLocator *vLocator=0 );
    6973  ~G4PropagatorInField();
    7074
    … …  
    116120   inline G4FieldTrack GetEndState() const;
    117121
    118    // The following methods are now obsolescent but *for now* will work
    119    //   They are being replaced by same-name methods in G4FieldManager,
    120    //   allowing the specialisation in different volumes.
    121    //   Their new behaviour is to change the values for the global field manager.
    122    inline G4double  GetMinimumEpsilonStep() const;
    123    inline void      SetMinimumEpsilonStep( G4double newEpsMin );
    124      // Minimum for Relative accuracy of any Step
    125 
     122   inline G4double  GetMinimumEpsilonStep() const;  // Min for relative accuracy
     123   inline void      SetMinimumEpsilonStep( G4double newEpsMin ); //  of any step
    126124   inline G4double  GetMaximumEpsilonStep() const;
    127125   inline void      SetMaximumEpsilonStep( G4double newEpsMax );
    128 
    129126   inline void      SetLargestAcceptableStep( G4double newBigDist );
    130127   inline G4double  GetLargestAcceptableStep();
     128     // The 6 above methods are now obsolescent but *for now* will work
     129     // They are being replaced by same-name methods in G4FieldManager,
     130     // allowing the specialisation in different volumes.
     131     // Their new behaviour is to change the values for the global field
     132     // manager
     133
     134   void SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter);
     135     // Set the filter that examines & stores 'intermediate'
     136     //  curved trajectory points.  Currently only position is stored.
     137
     138   std::vector<G4ThreeVector>* GimmeTrajectoryVectorAndForgetIt() const;
     139     // Access the points which have passed by the filter.
     140     // Responsibility for deleting the points lies with the client.
     141     // This method MUST BE called exactly ONCE per step.
     142
     143   void ClearPropagatorState();
     144     // Clear all the State of this class and its current associates
     145     //   --> the current field manager & chord finder will also be called
     146
     147   inline void SetDetectorFieldManager( G4FieldManager* newGlobalFieldManager );
     148     // Update this (dangerous) state -- for the time being
    131149 
    132      // Use alternative Locator(based on Brent Method,second order Intersection)
    133   inline void      SetBrentMethod(G4bool newLocator);
    134   inline G4bool    GetBrentMethod();
    135 
    136  public:  // with description
    137 
    138    // The following methods are obsolete and will not work --
    139    //   as they have been replaced by the same methods in G4FieldManager
    140    //   since Geant4 4.0
     150   inline void   SetUseSafetyForOptimization( G4bool );
     151   inline G4bool GetUseSafetyForOptimization();
     152     // Toggle & view parameter for using safety to discard
     153     //   unneccesary calls to navigator (thus 'optimising' performance)
     154   inline G4bool IntersectChord( G4ThreeVector  StartPointA,
     155                                 G4ThreeVector  EndPointB,
     156                                 G4double      &NewSafety,
     157                                 G4double      &LinearStepLength,
     158                                 G4ThreeVector &IntersectionPoint);
     159     // Intersect the chord from StartPointA to EndPointB
     160     // and return whether an intersection occurred
     161     // NOTE : SAFETY IS CHANGED
     162
     163   inline G4VIntersectionLocator* GetIntersectionLocator();
     164   inline void SetIntersectionLocator(G4VIntersectionLocator *pLocator );
     165 
     166 public:  // without description
     167
    141168   inline G4double  GetDeltaIntersection() const;
    142169   inline G4double  GetDeltaOneStep() const;
    … …  
    144171   inline void    SetDeltaIntersection( G4double deltaIntersection );
    145172   inline void    SetDeltaOneStep( G4double deltaOneStep ); 
    146 
    147  public:  // without description
     173     // The above 5 methods are obsolete and will not work, as they have been
     174     // replaced by the same methods in G4FieldManager since Geant4 4.0 ...
    148175
    149176   inline G4FieldManager*  GetCurrentFieldManager();
    150177   inline void             SetNavigatorForPropagating( G4Navigator *SimpleOrMultiNavigator );
    151178   inline G4Navigator*     GetNavigatorForPropagating();
    152 
    153  public:  // no description
    154179
    155180   inline void SetThresholdNoZeroStep( G4int noAct,
    … …  
    158183   inline G4int GetThresholdNoZeroSteps( G4int i );
    159184
    160  public:  // with description
    161   //
    162   void SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter);
    163   // Set the filter that examines & stores 'intermediate'
    164   //  curved trajectory points.  Currently only position is stored.
    165 
    166   std::vector<G4ThreeVector>* GimmeTrajectoryVectorAndForgetIt() const;
    167   // Access the points which have passed by the filter.
    168   // Responsibility for deleting the points lies with the client.
    169   // This method MUST BE called exactly ONCE per step.
    170 
    171   void ClearPropagatorState();
    172   // Clear all the State of this class and its current associates
    173   //   --> the current field manager & chord finder will also be called
    174 
    175   inline void SetDetectorFieldManager( G4FieldManager* newGlobalFieldManager );
    176       // Update this (dangerous) state -- for the time being
    177  
    178   inline void   SetUseSafetyForOptimization( G4bool );
    179   inline G4bool GetUseSafetyForOptimization();
    180       // Toggle & view parameter for using safety to discard
    181       //   unneccesary calls to navigator (thus 'optimising' performance)
    182 
    183185 protected:  // with description
    184 
    185    G4bool LocateIntersectionPoint(
    186         const  G4FieldTrack&       curveStartPointTangent,  //  A
    187         const  G4FieldTrack&       curveEndPointTangent,    //  B
    188         const  G4ThreeVector&      trialPoint,              //  E
    189                G4FieldTrack&       intersectPointTangent,   // Output
    190                G4bool&             recalculatedEndPoint);   // Out:
    191 
    192      // If such an intersection exists, this function
    193      // calculate the intersection point of the true path of the particle
    194      // with the surface of the current volume (or of one of its daughters).
    195      // (Should use lateral displacement as measure of convergence).
    196 
    197    G4bool IntersectChord( G4ThreeVector  StartPointA,
    198                           G4ThreeVector  EndPointB,
    199                           G4double      &NewSafety,
    200                           G4double      &LinearStepLength,
    201                           G4ThreeVector &IntersectionPoint);
    202      // Intersect the chord from StartPointA to EndPointB
    203      // and return whether an intersection occurred
    204 
    205    G4FieldTrack ReEstimateEndpoint( const G4FieldTrack &CurrentStateA, 
    206                                     const G4FieldTrack &EstimtdEndStateB,
    207                                           G4double      linearDistSq,
    208                                           G4double      curveDist);
    209      // Return new estimate for state after curveDist
    210      // starting from CurrentStateA,  to replace EstimtdEndStateB,
    211      // (and report displacement -- if field is compiled verbose.)
    212186
    213187   void PrintStepLengthDiagnostic( G4double      currentProposedStepLength,
    … …  
    217191 private:
    218192
    219   // ----------------------------------------------------------------------
    220   //  DATA Members
    221   // ----------------------------------------------------------------------
     193   // ----------------------------------------------------------------------
     194   //  DATA Members
     195   // ----------------------------------------------------------------------
    222196
    223197   G4FieldManager *fDetectorFieldMgr;
    … …  
    228202
    229203   G4Navigator   *fNavigator;
    230 
     204 
    231205   //  STATE information
    232206   //  -----------------
    … …  
    270244     // Geometrical tolerance defining surface thickness
    271245
    272   G4int maxNumberOfStepsForIntersection;
    273   G4int maxNumberOfCallsToReIntegration;
    274   G4int maxNumberOfCallsToReIntegration_depth;
    275     //  Counters for Statistics about Location and ReIntegrations
    276 private:
    277 
    278    static const G4int max_depth=4;
    279    G4FieldTrack* ptrInterMedFT[max_depth+1];
    280      // Used to store intermediate values of tracks in case of
    281      // too slow progress
    282 private:
    283    G4bool fUseBrentLocator;
    284    
    285 private:
    286 
    287   G4VCurvedTrajectoryFilter* fpTrajectoryFilter;
    288     // The filter encapsulates the algorithm which selects which
    289     // intermediate points should be stored in a trajectory.
    290     // When it is NULL, no intermediate points will be stored.
    291     // Else PIF::ComputeStep must submit (all) intermediate
    292     // points it calculates, to this filter.  (jacek 04/11/2002)
     246   G4VIntersectionLocator *fIntersectionLocator;
     247   G4bool fAllocatedLocator;
     248     // Used to Intersection Locator
     249
     250 private:
     251
     252   G4VCurvedTrajectoryFilter* fpTrajectoryFilter;
     253     // The filter encapsulates the algorithm which selects which
     254     // intermediate points should be stored in a trajectory.
     255     // When it is NULL, no intermediate points will be stored.
     256     // Else PIF::ComputeStep must submit (all) intermediate
     257     // points it calculates, to this filter.  (jacek 04/11/2002)
    293258};
    294259

  • trunk/source/geometry/navigation/include/G4PropagatorInField.icc

    r850 r921  
    2525//
    2626//
    27 // $Id: G4PropagatorInField.icc,v 1.11 2008/05/28 09:12:05 tnikitin Exp $
    28 // GEANT4 tag $Name: HEAD $
     27// $Id: G4PropagatorInField.icc,v 1.13 2008/10/29 14:31:55 gcosmo Exp $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//
    … …  
    268268}
    269269
    270 inline
    271 void G4PropagatorInField::SetBrentMethod(G4bool newLocator)
    272 {
    273   fUseBrentLocator=newLocator;
    274 }
    275 inline
    276 G4bool G4PropagatorInField::GetBrentMethod()
    277 {
    278   return fUseBrentLocator;
    279 }
     270inline
     271void G4PropagatorInField::
     272SetIntersectionLocator( G4VIntersectionLocator *pIntLoc )
     273{
     274  if(pIntLoc)  { fIntersectionLocator= pIntLoc; }
     275}
     276
     277inline
     278G4VIntersectionLocator* G4PropagatorInField::GetIntersectionLocator()
     279{
     280  return fIntersectionLocator;
     281}
     282
     283inline
     284G4bool G4PropagatorInField::IntersectChord( G4ThreeVector  StartPointA,
     285                                            G4ThreeVector  EndPointB,
     286                                            G4double      &NewSafety,
     287                                            G4double      &LinearStepLength,
     288                                            G4ThreeVector &IntersectionPoint )
     289{
     290  // Calculate the direction and length of the chord AB
     291  //
     292  return fIntersectionLocator
     293         ->IntersectChord(StartPointA,EndPointB,NewSafety,
     294                          fPreviousSafety,fPreviousSftOrigin,
     295                          LinearStepLength,IntersectionPoint);
     296}

  • trunk/source/geometry/navigation/include/G4RegularNavigation.hh

    r850 r921  
    2626//
    2727// $Id: G4RegularNavigation.hh,v 1.2 2007/10/18 14:18:36 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4ReplicaNavigation.hh

    r850 r921  
    2626//
    2727// $Id: G4ReplicaNavigation.hh,v 1.6 2007/05/18 07:31:09 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4ReplicaNavigation.icc

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

  • trunk/source/geometry/navigation/include/G4SafetyHelper.hh

    r850 r921  
    2626//
    2727// $Id: G4SafetyHelper.hh,v 1.7 2007/05/02 15:32:13 japost Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030//

  • trunk/source/geometry/navigation/include/G4TransportationManager.hh

    r850 r921  
    2626//
    2727// $Id: G4TransportationManager.hh,v 1.12 2007/04/20 15:28:37 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929//
    3030// class G4TransportationManager

  • trunk/source/geometry/navigation/include/G4TransportationManager.icc

    r850 r921  
    2626//
    2727// $Id: G4TransportationManager.icc,v 1.10 2007/04/20 15:28:37 gcosmo Exp $
    28 // GEANT4 tag $Name: HEAD $
     28// GEANT4 tag $Name: geant4-09-02-cand-01 $
    2929// ------------------------------------------------------------
    3030//  GEANT 4  inlined function members implementation

  • trunk/source/geometry/navigation/include/G4VoxelNavigation.hh

    r850 r921  
    2626//
    2727// $Id: G4VoxelNavigation.hh,v 1.5 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/include/G4VoxelNavigation.icc

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

  • 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.