| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4SimpleLocator.cc 90451 2015-05-29 09:48:07Z gcosmo $
>> 27 //
26 // Class G4SimpleLocator implementation 28 // Class G4SimpleLocator implementation
27 // 29 //
28 // 27.10.08 - Tatiana Nikitina, extracted from 30 // 27.10.08 - Tatiana Nikitina, extracted from G4PropagatorInField class
29 // 04.10.11 - John Apostolakis, revised conver 31 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
30 // ------------------------------------------- 32 // ---------------------------------------------------------------------------
31 33
32 #include <iomanip> 34 #include <iomanip>
33 35
34 #include "G4ios.hh" 36 #include "G4ios.hh"
35 #include "G4SimpleLocator.hh" 37 #include "G4SimpleLocator.hh"
36 38
37 G4SimpleLocator::G4SimpleLocator(G4Navigator * 39 G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator)
38 : G4VIntersectionLocator(theNavigator) 40 : G4VIntersectionLocator(theNavigator)
39 { 41 {
40 } 42 }
41 43
42 G4SimpleLocator::~G4SimpleLocator() = default; << 44 G4SimpleLocator::~G4SimpleLocator()
>> 45 {
>> 46 }
43 47
44 // ------------------------------------------- 48 // --------------------------------------------------------------------------
45 // G4bool G4PropagatorInField::LocateIntersect 49 // G4bool G4PropagatorInField::LocateIntersectionPoint(
46 // const G4FieldTrack& CurveStart 50 // const G4FieldTrack& CurveStartPointVelocity, // A
47 // const G4FieldTrack& CurveEndPo 51 // const G4FieldTrack& CurveEndPointVelocity, // B
48 // const G4ThreeVector& TrialPoint 52 // const G4ThreeVector& TrialPoint, // E
49 // G4FieldTrack& Intersecte 53 // G4FieldTrack& IntersectedOrRecalculated // Output
50 // G4bool& recalculat 54 // G4bool& recalculated ) // Out
51 // ------------------------------------------- 55 // --------------------------------------------------------------------------
52 // 56 //
53 // Function that returns the intersection of t 57 // Function that returns the intersection of the true path with the surface
54 // of the current volume (either the external 58 // of the current volume (either the external one or the inner one with one
55 // of the daughters: 59 // of the daughters:
56 // 60 //
57 // A = Initial point 61 // A = Initial point
58 // B = another point 62 // B = another point
59 // 63 //
60 // Both A and B are assumed to be on the true 64 // Both A and B are assumed to be on the true path:
61 // 65 //
62 // E is the first point of intersection of 66 // E is the first point of intersection of the chord AB with
63 // a volume other than A (on the surface 67 // a volume other than A (on the surface of A or of a daughter)
64 // 68 //
65 // Convention of Use : 69 // Convention of Use :
66 // i) If it returns "true", then Intersect 70 // i) If it returns "true", then IntersectionPointVelocity is set
67 // to the approximate intersection point 71 // to the approximate intersection point.
68 // ii) If it returns "false", no intersecti 72 // ii) If it returns "false", no intersection was found.
69 // The validity of IntersectedOrRecal 73 // The validity of IntersectedOrRecalculated depends on 'recalculated'
70 // a) if latter is false, then Intersec 74 // a) if latter is false, then IntersectedOrRecalculated is invalid.
71 // b) if latter is true, then Intersec 75 // b) if latter is true, then IntersectedOrRecalculated is
72 // the new endpoint, due to a re-i 76 // the new endpoint, due to a re-integration.
73 // ------------------------------------------- 77 // --------------------------------------------------------------------------
74 // NOTE: implementation taken from G4Propagato 78 // NOTE: implementation taken from G4PropagatorInField
75 // 79 //
76 G4bool G4SimpleLocator::EstimateIntersectionPo 80 G4bool G4SimpleLocator::EstimateIntersectionPoint(
77 const G4FieldTrack& CurveStart 81 const G4FieldTrack& CurveStartPointVelocity, // A
78 const G4FieldTrack& CurveEndPo 82 const G4FieldTrack& CurveEndPointVelocity, // B
79 const G4ThreeVector& TrialPoint 83 const G4ThreeVector& TrialPoint, // E
80 G4FieldTrack& Intersecte 84 G4FieldTrack& IntersectedOrRecalculatedFT, // Output
81 G4bool& recalculat 85 G4bool& recalculatedEndPoint, // Out
82 G4double &fPrevious 86 G4double &fPreviousSafety, //In/Out
83 G4ThreeVector &fPrevious 87 G4ThreeVector &fPreviousSftOrigin) //In/Out
84 { 88 {
85 // Find Intersection Point ( A, B, E ) of t 89 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
86 90
87 G4bool found_approximate_intersection = fals 91 G4bool found_approximate_intersection = false;
88 G4bool there_is_no_intersection = fals 92 G4bool there_is_no_intersection = false;
89 93
90 G4FieldTrack CurrentA_PointVelocity = Curve 94 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
91 G4FieldTrack CurrentB_PointVelocity = Curve 95 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
92 G4ThreeVector CurrentE_Point = TrialPoint; 96 G4ThreeVector CurrentE_Point = TrialPoint;
93 G4bool validNormalAtE = false; 97 G4bool validNormalAtE = false;
94 G4ThreeVector NormalAtEntry; 98 G4ThreeVector NormalAtEntry;
95 99
96 G4FieldTrack ApproxIntersecPointV(CurveEndP 100 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
97 G4double NewSafety = 0.0; 101 G4double NewSafety = 0.0;
98 G4bool last_AF_intersection = false; 102 G4bool last_AF_intersection = false;
99 G4bool final_section = true; // Shows wheth 103 G4bool final_section = true; // Shows whether current section is last
100 // (i.e. B=ful 104 // (i.e. B=full end)
101 recalculatedEndPoint = false; 105 recalculatedEndPoint = false;
102 106
103 G4bool restoredFullEndpoint = false; 107 G4bool restoredFullEndpoint = false;
104 108
105 G4int substep_no = 0; 109 G4int substep_no = 0;
106 110
107 // Limits for substep number 111 // Limits for substep number
108 // 112 //
109 const G4int max_substeps = 100000000; // T 113 const G4int max_substeps = 100000000; // Test 120 (old value 100 )
110 const G4int warn_substeps = 1000; // 114 const G4int warn_substeps = 1000; // 100
111 115
112 // Statistics for substeps 116 // Statistics for substeps
113 // 117 //
114 static G4ThreadLocal G4int max_no_seen= -1; 118 static G4ThreadLocal G4int max_no_seen= -1;
115 119
116 NormalAtEntry = GetSurfaceNormal( CurrentE_P 120 NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE);
117 121
118 #ifdef G4DEBUG_FIELD 122 #ifdef G4DEBUG_FIELD
119 const G4double tolerance = 1.0e-8; << 123 static G4double tolerance = 1.0e-8;
120 G4ThreeVector StartPosition= CurveStartPoin 124 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
121 if( (TrialPoint - StartPosition).mag() < tol 125 if( (TrialPoint - StartPosition).mag() < tolerance * CLHEP::mm )
122 { 126 {
123 G4Exception("G4SimpleLocator::EstimateInt 127 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
124 "GeomNav1002", JustWarning, 128 "GeomNav1002", JustWarning,
125 "Intersection point F is exac 129 "Intersection point F is exactly at start point A." );
126 } 130 }
127 #endif 131 #endif
128 132
129 do 133 do
130 { 134 {
131 G4ThreeVector Point_A = CurrentA_PointVel 135 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
132 G4ThreeVector Point_B = CurrentB_PointVel 136 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
133 137
134 // F = a point on true AB path close to p 138 // F = a point on true AB path close to point E
135 // (the closest if possible) 139 // (the closest if possible)
136 // 140 //
137 ApproxIntersecPointV = GetChordFinderFor( 141 ApproxIntersecPointV = GetChordFinderFor()
138 ->ApproxCurvePoint 142 ->ApproxCurvePointV( CurrentA_PointVelocity,
139 143 CurrentB_PointVelocity,
140 144 CurrentE_Point,
141 145 GetEpsilonStepFor());
142 // The above method is the key & most 146 // The above method is the key & most intuitive part ...
143 147
144 #ifdef G4DEBUG_FIELD 148 #ifdef G4DEBUG_FIELD
145 if( ApproxIntersecPointV.GetCurveLength() 149 if( ApproxIntersecPointV.GetCurveLength() >
146 CurrentB_PointVelocity.GetCurveLength 150 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
147 { 151 {
148 G4Exception("G4SimpleLocator::EstimateI 152 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
149 "GeomNav0003", FatalExcepti 153 "GeomNav0003", FatalException,
150 "Intermediate F point is pa 154 "Intermediate F point is past end B point!" );
151 } 155 }
152 #endif 156 #endif
153 157
154 G4ThreeVector CurrentF_Point = ApproxInte << 158 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
155 159
156 // First check whether EF is small - then 160 // First check whether EF is small - then F is a good approx. point
157 // Calculate the length and direction of 161 // Calculate the length and direction of the chord AF
158 // 162 //
159 G4ThreeVector ChordEF_Vector = CurrentF_P << 163 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
160 164
161 G4ThreeVector NewMomentumDir = ApproxInte << 165 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
162 G4double MomDir_dot_Norm = NewMomentumDir << 166 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
163 167
164 G4ThreeVector ChordAB = Point_B - Point_A << 168 G4ThreeVector ChordAB = Point_B - Point_A;
165 169
166 #ifdef G4DEBUG_FIELD 170 #ifdef G4DEBUG_FIELD
167 G4VIntersectionLocator:: 171 G4VIntersectionLocator::
168 ReportTrialStep( substep_no, ChordAB, C 172 ReportTrialStep( substep_no, ChordAB, ChordEF_Vector,
169 NewMomentumDir, Normal << 173 NewMomentumDir, NormalAtEntry, validNormalAtE );
170 #endif 174 #endif
171 // Check Sign is always exiting !! TODO 175 // Check Sign is always exiting !! TODO
172 // Could ( > -epsilon) be used instead? 176 // Could ( > -epsilon) be used instead?
173 // 177 //
174 G4bool adequate_angle = ( MomDir_dot_Norm 178 G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 )
175 || (! validNormalAtE 179 || (! validNormalAtE ); // Invalid
176 G4double EF_dist2= ChordEF_Vector.mag2(); 180 G4double EF_dist2= ChordEF_Vector.mag2();
177 if ( ( EF_dist2 <= sqr(fiDeltaIntersecti 181 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
178 || ( EF_dist2 <= kCarTolerance*kCarTole 182 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
179 { 183 {
180 found_approximate_intersection = true; 184 found_approximate_intersection = true;
181 185
182 // Create the "point" return value 186 // Create the "point" return value
183 // 187 //
184 IntersectedOrRecalculatedFT = ApproxInt 188 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
185 IntersectedOrRecalculatedFT.SetPosition 189 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
186 190
187 if ( GetAdjustementOfFoundIntersection( 191 if ( GetAdjustementOfFoundIntersection() )
188 { 192 {
189 // Try to Get Correction of Intersect 193 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
190 // 194 //
191 G4ThreeVector IP; 195 G4ThreeVector IP;
192 G4ThreeVector MomentumDir= ApproxInte 196 G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
193 G4bool goodCorrection = AdjustmentOfF 197 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
194 CurrentE_Po 198 CurrentE_Point, CurrentF_Point, MomentumDir,
195 last_AF_int 199 last_AF_intersection, IP, NewSafety,
196 fPreviousSa 200 fPreviousSafety, fPreviousSftOrigin );
197 201
198 if ( goodCorrection ) << 202 if(goodCorrection)
199 { 203 {
200 IntersectedOrRecalculatedFT = Appro 204 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
201 IntersectedOrRecalculatedFT.SetPosi 205 IntersectedOrRecalculatedFT.SetPosition(IP);
202 } 206 }
203 } 207 }
204 208
205 // Note: in order to return a point on 209 // Note: in order to return a point on the boundary,
206 // we must return E. But it is F 210 // we must return E. But it is F on the curve.
207 // So we must "cheat": we are usi 211 // So we must "cheat": we are using the position at point E
208 // and the velocity at point F !! 212 // and the velocity at point F !!!
209 // 213 //
210 // This must limit the length we can al 214 // This must limit the length we can allow for displacement!
211 } 215 }
212 else // E is NOT close enough to the cur 216 else // E is NOT close enough to the curve (ie point F)
213 { 217 {
214 // Check whether any volumes are encoun 218 // Check whether any volumes are encountered by the chord AF
215 // ------------------------------------ 219 // ---------------------------------------------------------
216 // First relocate to restore any Voxel 220 // First relocate to restore any Voxel etc information
217 // in the Navigator before calling Comp 221 // in the Navigator before calling ComputeStep()
218 // 222 //
219 GetNavigatorFor()->LocateGlobalPointWit 223 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
220 224
221 G4ThreeVector PointG; // Candidate in 225 G4ThreeVector PointG; // Candidate intersection point
222 G4double stepLengthAF; 226 G4double stepLengthAF;
223 G4bool usedNavigatorAF = false; 227 G4bool usedNavigatorAF = false;
224 G4bool Intersects_AF = IntersectChord( 228 G4bool Intersects_AF = IntersectChord( Point_A,
225 229 CurrentF_Point,
226 230 NewSafety,
227 231 fPreviousSafety,
228 232 fPreviousSftOrigin,
229 233 stepLengthAF,
230 234 PointG,
231 235 &usedNavigatorAF );
232 last_AF_intersection = Intersects_AF; 236 last_AF_intersection = Intersects_AF;
233 if( Intersects_AF ) 237 if( Intersects_AF )
234 { 238 {
235 // G is our new Candidate for the int 239 // G is our new Candidate for the intersection point.
236 // It replaces "E" and we will repea 240 // It replaces "E" and we will repeat the test to see if
237 // it is a good enough approximate po 241 // it is a good enough approximate point for us.
238 // B <- F 242 // B <- F
239 // E <- G 243 // E <- G
240 244
241 CurrentB_PointVelocity = ApproxInters 245 CurrentB_PointVelocity = ApproxIntersecPointV;
242 CurrentE_Point = PointG; 246 CurrentE_Point = PointG;
243 247
244 // Need to recalculate the Exit Norma 248 // Need to recalculate the Exit Normal at the new PointG
245 // Relies on a call to Navigator::Com 249 // Relies on a call to Navigator::ComputeStep in IntersectChord above
246 // If the safety was adequate (for th 250 // If the safety was adequate (for the step) this would NOT be called!
247 // But this must not occur, no inters 251 // But this must not occur, no intersection can be found in that case,
248 // so this branch, ie if( Intersects_ 252 // so this branch, ie if( Intersects_AF ) would not be reached!
249 // 253 //
250 G4bool validNormalLast; 254 G4bool validNormalLast;
251 NormalAtEntry = GetSurfaceNormal( Po 255 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
252 validNormalAtE = validNormalLast; 256 validNormalAtE = validNormalLast;
253 257
254 // By moving point B, must take care 258 // By moving point B, must take care if current
255 // AF has no intersection to try curr 259 // AF has no intersection to try current FB!!
256 // 260 //
257 final_section = false; << 261 final_section= false;
258 262
259 #ifdef G4VERBOSE 263 #ifdef G4VERBOSE
260 if( fVerboseLevel > 3 ) 264 if( fVerboseLevel > 3 )
261 { 265 {
262 G4cout << "G4PiF::LI> Investigating 266 G4cout << "G4PiF::LI> Investigating intermediate point"
263 << " at s=" << ApproxInterse 267 << " at s=" << ApproxIntersecPointV.GetCurveLength()
264 << " on way to full s=" 268 << " on way to full s="
265 << CurveEndPointVelocity.Get 269 << CurveEndPointVelocity.GetCurveLength() << G4endl;
266 } 270 }
267 #endif 271 #endif
268 } 272 }
269 else // not Intersects_AF 273 else // not Intersects_AF
270 { 274 {
271 // In this case: 275 // In this case:
272 // There is NO intersection of AF wit 276 // There is NO intersection of AF with a volume boundary.
273 // We must continue the search in the 277 // We must continue the search in the segment FB!
274 // 278 //
275 GetNavigatorFor()->LocateGlobalPointW 279 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
276 280
277 G4double stepLengthFB; 281 G4double stepLengthFB;
278 G4ThreeVector PointH; 282 G4ThreeVector PointH;
279 G4bool usedNavigatorFB = false; << 283 G4bool usedNavigatorFB=false;
280 284
281 // Check whether any volumes are enco 285 // Check whether any volumes are encountered by the chord FB
282 // ---------------------------------- 286 // ---------------------------------------------------------
283 287
284 G4bool Intersects_FB = IntersectChord 288 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
285 289 NewSafety,fPreviousSafety,
286 290 fPreviousSftOrigin,
287 291 stepLengthFB,
288 292 PointH, &usedNavigatorFB );
289 if( Intersects_FB ) 293 if( Intersects_FB )
290 { 294 {
291 // There is an intersection of FB w 295 // There is an intersection of FB with a volume boundary
292 // H <- First Intersection of Chord 296 // H <- First Intersection of Chord FB
293 297
294 // H is our new Candidate for the i 298 // H is our new Candidate for the intersection point.
295 // It replaces "E" and we will rep 299 // It replaces "E" and we will repeat the test to see if
296 // it is a good enough approximate 300 // it is a good enough approximate point for us.
297 301
298 // Note that F must be in volume vo 302 // Note that F must be in volume volA (the same as A)
299 // (otherwise AF would meet a volum 303 // (otherwise AF would meet a volume boundary!)
300 // A <- F 304 // A <- F
301 // E <- H 305 // E <- H
302 // 306 //
303 CurrentA_PointVelocity = ApproxInte 307 CurrentA_PointVelocity = ApproxIntersecPointV;
304 CurrentE_Point = PointH; 308 CurrentE_Point = PointH;
305 309
306 // Need to recalculate the Exit Nor 310 // Need to recalculate the Exit Normal at the new PointG
307 // Relies on call to Navigator::Com 311 // Relies on call to Navigator::ComputeStep in IntersectChord above
308 // If safety was adequate (for the 312 // If safety was adequate (for the step) this would NOT be called!
309 // But this must not occur, no inte 313 // But this must not occur, no intersection found in that case,
310 // so this branch, i.e. if( Interse 314 // so this branch, i.e. if( Intersects_AF ) would not be reached!
311 // 315 //
312 G4bool validNormalLast; 316 G4bool validNormalLast;
313 NormalAtEntry = GetSurfaceNormal( 317 NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast );
314 validNormalAtE = validNormalLast; 318 validNormalAtE = validNormalLast;
315 } 319 }
316 else // not Intersects_FB 320 else // not Intersects_FB
317 { 321 {
318 // There is NO intersection of FB w 322 // There is NO intersection of FB with a volume boundary
319 323
320 if( final_section ) 324 if( final_section )
321 { 325 {
322 // If B is the original endpoint, 326 // If B is the original endpoint, this means that whatever
323 // volume(s) intersected the orig 327 // volume(s) intersected the original chord, none touch the
324 // smaller chords we have used. 328 // smaller chords we have used.
325 // The value of 'IntersectedOrRec 329 // The value of 'IntersectedOrRecalculatedFT' returned is
326 // likely not valid 330 // likely not valid
327 331
328 there_is_no_intersection = true; 332 there_is_no_intersection = true; // real final_section
329 } 333 }
330 else 334 else
331 { 335 {
332 // We must restore the original e 336 // We must restore the original endpoint
333 337
334 CurrentA_PointVelocity = CurrentB 338 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
335 CurrentB_PointVelocity = CurveEnd 339 CurrentB_PointVelocity = CurveEndPointVelocity;
336 restoredFullEndpoint = true; 340 restoredFullEndpoint = true;
337 } 341 }
338 } // Endif (Intersects_FB) 342 } // Endif (Intersects_FB)
339 } // Endif (Intersects_AF) 343 } // Endif (Intersects_AF)
340 344
341 // Ensure that the new endpoints are no 345 // Ensure that the new endpoints are not further apart in space
342 // than on the curve due to different e 346 // than on the curve due to different errors in the integration
343 // 347 //
344 G4double linDistSq, curveDist; 348 G4double linDistSq, curveDist;
345 linDistSq = ( CurrentB_PointVelocity.Ge 349 linDistSq = ( CurrentB_PointVelocity.GetPosition()
346 - CurrentA_PointVelocity.Ge 350 - CurrentA_PointVelocity.GetPosition() ).mag2();
347 curveDist = CurrentB_PointVelocity.GetC 351 curveDist = CurrentB_PointVelocity.GetCurveLength()
348 - CurrentA_PointVelocity.Ge 352 - CurrentA_PointVelocity.GetCurveLength();
349 353
350 // Change this condition for very stric 354 // Change this condition for very strict parameters of propagation
351 // 355 //
352 if( curveDist*curveDist*(1+2* GetEpsilo 356 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
353 { 357 {
354 // Re-integrate to obtain a new B 358 // Re-integrate to obtain a new B
355 // 359 //
356 G4FieldTrack newEndPointFT = 360 G4FieldTrack newEndPointFT =
357 ReEstimateEndpoint( CurrentA_ 361 ReEstimateEndpoint( CurrentA_PointVelocity,
358 CurrentB_ 362 CurrentB_PointVelocity,
359 linDistSq 363 linDistSq, // to avoid recalculation
360 curveDist 364 curveDist );
361 G4FieldTrack oldPointVelB = CurrentB_ 365 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
362 CurrentB_PointVelocity = newEndPointF 366 CurrentB_PointVelocity = newEndPointFT;
363 367
364 if( (final_section)) // real final se 368 if( (final_section)) // real final section
365 { 369 {
366 recalculatedEndPoint = true; 370 recalculatedEndPoint = true;
367 IntersectedOrRecalculatedFT = newEn 371 IntersectedOrRecalculatedFT = newEndPointFT;
368 // So that we can return it, if i 372 // So that we can return it, if it is the endpoint!
369 } 373 }
370 } 374 }
371 if( curveDist < 0.0 ) 375 if( curveDist < 0.0 )
372 { 376 {
373 fVerboseLevel = 5; // Print out a max 377 fVerboseLevel = 5; // Print out a maximum of information
374 printStatus( CurrentA_PointVelocity, 378 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
375 -1.0, NewSafety, subste 379 -1.0, NewSafety, substep_no );
376 std::ostringstream message; 380 std::ostringstream message;
377 message << "Error in advancing propag 381 message << "Error in advancing propagation." << G4endl
378 << " Point A (start) i 382 << " Point A (start) is " << CurrentA_PointVelocity
379 << G4endl 383 << G4endl
380 << " Point B (end) i 384 << " Point B (end) is " << CurrentB_PointVelocity
381 << G4endl 385 << G4endl
382 << " Curve distance is 386 << " Curve distance is " << curveDist << G4endl
383 << G4endl 387 << G4endl
384 << "The final curve point is 388 << "The final curve point is not further along"
385 << " than the original!" << G 389 << " than the original!" << G4endl;
386 390
387 if( recalculatedEndPoint ) 391 if( recalculatedEndPoint )
388 { 392 {
389 message << "Recalculation of EndPoi 393 message << "Recalculation of EndPoint was called with fEpsStep= "
390 << GetEpsilonStepFor() << G 394 << GetEpsilonStepFor() << G4endl;
391 } 395 }
392 message.precision(20); 396 message.precision(20);
393 message << " Point A (Curve start) 397 message << " Point A (Curve start) is " << CurveStartPointVelocity
394 << G4endl 398 << G4endl
395 << " Point B (Curve end) 399 << " Point B (Curve end) is " << CurveEndPointVelocity
396 << G4endl 400 << G4endl
397 << " Point A (Current start) 401 << " Point A (Current start) is " << CurrentA_PointVelocity
398 << G4endl 402 << G4endl
399 << " Point B (Current end) 403 << " Point B (Current end) is " << CurrentB_PointVelocity
400 << G4endl 404 << G4endl
401 << " Point E (Trial Point) 405 << " Point E (Trial Point) is " << CurrentE_Point
402 << G4endl 406 << G4endl
403 << " Point F (Intersection) 407 << " Point F (Intersection) is " << ApproxIntersecPointV
404 << G4endl 408 << G4endl
405 << " LocateIntersectio 409 << " LocateIntersection parameters are : Substep no= "
406 << substep_no; 410 << substep_no;
407 411
408 G4Exception("G4SimpleLocator::Estimat 412 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
409 "GeomNav0003", FatalExcep 413 "GeomNav0003", FatalException, message);
410 } 414 }
411 415
412 if ( restoredFullEndpoint ) << 416 if(restoredFullEndpoint)
413 { 417 {
414 final_section = restoredFullEndpoint; 418 final_section = restoredFullEndpoint;
415 restoredFullEndpoint = false; 419 restoredFullEndpoint = false;
416 } 420 }
417 } // EndIf ( E is close enough to the cur 421 } // EndIf ( E is close enough to the curve, ie point F. )
418 // tests ChordAF_Vector.mag() <= maximu 422 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
419 423
420 #ifdef G4DEBUG_LOCATE_INTERSECTION 424 #ifdef G4DEBUG_LOCATE_INTERSECTION
421 G4int trigger_substepno_print= warn_subst << 425 static G4int trigger_substepno_print= warn_substeps - 20;
422 426
423 if( substep_no >= trigger_substepno_print 427 if( substep_no >= trigger_substepno_print )
424 { 428 {
425 G4cout << "Difficulty in converging in 429 G4cout << "Difficulty in converging in "
426 << "G4SimpleLocator::EstimateInt 430 << "G4SimpleLocator::EstimateIntersectionPoint():"
427 << G4endl 431 << G4endl
428 << " Substep no = " << subste 432 << " Substep no = " << substep_no << G4endl;
429 if( substep_no == trigger_substepno_pri 433 if( substep_no == trigger_substepno_print )
430 { 434 {
431 printStatus( CurveStartPointVelocity, 435 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
432 -1.0, NewSafety, 0); 436 -1.0, NewSafety, 0);
433 } 437 }
434 G4cout << " State of point A: "; 438 G4cout << " State of point A: ";
435 printStatus( CurrentA_PointVelocity, Cu 439 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
436 -1.0, NewSafety, substep_n 440 -1.0, NewSafety, substep_no-1, 0);
437 G4cout << " State of point B: "; 441 G4cout << " State of point B: ";
438 printStatus( CurrentA_PointVelocity, Cu 442 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
439 -1.0, NewSafety, substep_n 443 -1.0, NewSafety, substep_no);
440 } 444 }
441 #endif 445 #endif
442 ++substep_no; << 446 substep_no++;
443 447
444 } while ( ( ! found_approximate_intersection 448 } while ( ( ! found_approximate_intersection )
445 && ( ! there_is_no_intersection ) 449 && ( ! there_is_no_intersection )
446 && ( substep_no <= max_substeps) ); 450 && ( substep_no <= max_substeps) ); // UNTIL found or failed
447 451
448 if( substep_no > max_no_seen ) 452 if( substep_no > max_no_seen )
449 { 453 {
450 max_no_seen = substep_no; 454 max_no_seen = substep_no;
451 #ifdef G4DEBUG_LOCATE_INTERSECTION 455 #ifdef G4DEBUG_LOCATE_INTERSECTION
452 if( max_no_seen > warn_substeps ) 456 if( max_no_seen > warn_substeps )
453 { 457 {
454 trigger_substepno_print = max_no_seen-20 458 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
455 } 459 }
456 #endif 460 #endif
457 } 461 }
458 462
459 if( ( substep_no >= max_substeps) 463 if( ( substep_no >= max_substeps)
460 && !there_is_no_intersection 464 && !there_is_no_intersection
461 && !found_approximate_intersection ) 465 && !found_approximate_intersection )
462 { 466 {
463 G4cout << "ERROR - G4SimpleLocator::Estima 467 G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl
464 << " Start and Endpoint of R 468 << " Start and Endpoint of Requested Step:" << G4endl;
465 printStatus( CurveStartPointVelocity, Curv 469 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
466 -1.0, NewSafety, 0); 470 -1.0, NewSafety, 0);
467 G4cout << G4endl 471 G4cout << G4endl
468 << " Start and end-point of 472 << " Start and end-point of current Sub-Step:" << G4endl;
469 printStatus( CurrentA_PointVelocity, Curre 473 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
470 -1.0, NewSafety, substep_no-1 474 -1.0, NewSafety, substep_no-1);
471 printStatus( CurrentA_PointVelocity, Curre 475 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
472 -1.0, NewSafety, substep_no); 476 -1.0, NewSafety, substep_no);
473 477
474 std::ostringstream message; 478 std::ostringstream message;
475 message << "Convergence is requiring too m 479 message << "Convergence is requiring too many substeps: "
476 << substep_no << G4endl 480 << substep_no << G4endl
477 << " Abandoning effort to 481 << " Abandoning effort to intersect." << G4endl
478 << " Found intersection = 482 << " Found intersection = "
479 << found_approximate_intersection 483 << found_approximate_intersection << G4endl
480 << " Intersection exists 484 << " Intersection exists = "
481 << !there_is_no_intersection << G4 485 << !there_is_no_intersection << G4endl;
482 message.precision(10); 486 message.precision(10);
483 G4double done_len = CurrentA_PointVelocity 487 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
484 G4double full_len = CurveEndPointVelocity. 488 G4double full_len = CurveEndPointVelocity.GetCurveLength();
485 message << " Undertaken only leng 489 message << " Undertaken only length: " << done_len
486 << " out of " << full_len << " req 490 << " out of " << full_len << " required." << G4endl
487 << " Remaining length = " 491 << " Remaining length = " << full_len-done_len;
488 492
489 G4Exception("G4SimpleLocator::EstimateInte 493 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
490 "GeomNav0003", FatalException, 494 "GeomNav0003", FatalException, message);
491 } 495 }
492 else if( substep_no >= warn_substeps ) 496 else if( substep_no >= warn_substeps )
493 { 497 {
494 std::ostringstream message; 498 std::ostringstream message;
495 message.precision(10); 499 message.precision(10);
496 500
497 message << "Many substeps while trying to 501 message << "Many substeps while trying to locate intersection." << G4endl
498 << " Undertaken length: " 502 << " Undertaken length: "
499 << CurrentB_PointVelocity.GetCurve 503 << CurrentB_PointVelocity.GetCurveLength()
500 << " - Needed: " << substep_no << 504 << " - Needed: " << substep_no << " substeps." << G4endl
501 << " Warning level = " << 505 << " Warning level = " << warn_substeps
502 << " and maximum substeps = " << m 506 << " and maximum substeps = " << max_substeps;
503 G4Exception("G4SimpleLocator::EstimateInte 507 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
504 "GeomNav1002", JustWarning, me 508 "GeomNav1002", JustWarning, message);
505 } 509 }
506 return !there_is_no_intersection; // Succe 510 return !there_is_no_intersection; // Success or failure
507 } 511 }
508 512