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