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