| 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: G4MultiLevelLocator.cc 66872 2013-01-15 01:25:57Z japost $
>> 27 //
26 // Class G4MultiLevelLocator implementation 28 // Class G4MultiLevelLocator implementation
27 // 29 //
28 // 27.10.08 - Tatiana Nikitina. 30 // 27.10.08 - Tatiana Nikitina.
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 "G4MultiLevelLocator.hh" 37 #include "G4MultiLevelLocator.hh"
36 #include "G4LocatorChangeRecord.hh" <<
37 #include "G4LocatorChangeLogger.hh" <<
38 38
39 G4MultiLevelLocator::G4MultiLevelLocator(G4Nav 39 G4MultiLevelLocator::G4MultiLevelLocator(G4Navigator *theNavigator)
40 : G4VIntersectionLocator(theNavigator) << 40 : G4VIntersectionLocator(theNavigator)
41 { 41 {
42 // In case of too slow progress in finding I 42 // In case of too slow progress in finding Intersection Point
43 // intermediates Points on the Track must be 43 // intermediates Points on the Track must be stored.
44 // Initialise the array of Pointers [max_dep 44 // Initialise the array of Pointers [max_depth+1] to do this
45 45
46 G4ThreeVector zeroV(0.0,0.0,0.0); 46 G4ThreeVector zeroV(0.0,0.0,0.0);
47 for (auto & idepth : ptrInterMedFT) << 47 for (G4int idepth=0; idepth<max_depth+1; idepth++ )
48 { 48 {
49 idepth = new G4FieldTrack( zeroV, zeroV, 0 << 49 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
50 } 50 }
51 <<
52 if (fCheckMode) <<
53 { <<
54 // Trial values Loose Mediu <<
55 // To happen: Infrequent Occas <<
56 SetMaxSteps(150); // 300 85 <<
57 SetWarnSteps(80); // 250 60 <<
58 } <<
59 } 51 }
60 52
61 G4MultiLevelLocator::~G4MultiLevelLocator() 53 G4MultiLevelLocator::~G4MultiLevelLocator()
62 { 54 {
63 for (auto & idepth : ptrInterMedFT) << 55 for ( G4int idepth=0; idepth<max_depth+1; idepth++)
64 { 56 {
65 delete idepth; << 57 delete ptrInterMedFT[idepth];
66 } 58 }
67 #ifdef G4DEBUG_FIELD <<
68 ReportStatistics(); <<
69 #endif <<
70 } 59 }
71 60
72 <<
73 // ------------------------------------------- 61 // --------------------------------------------------------------------------
74 // G4bool G4PropagatorInField::LocateIntersect 62 // G4bool G4PropagatorInField::LocateIntersectionPoint(
75 // const G4FieldTrack& CurveStart 63 // const G4FieldTrack& CurveStartPointVelocity, // A
76 // const G4FieldTrack& CurveEndPo 64 // const G4FieldTrack& CurveEndPointVelocity, // B
77 // const G4ThreeVector& TrialPoint 65 // const G4ThreeVector& TrialPoint, // E
78 // G4FieldTrack& Intersecte 66 // G4FieldTrack& IntersectedOrRecalculated // Output
79 // G4bool& recalculat 67 // G4bool& recalculated ) // Out
80 // ------------------------------------------- 68 // --------------------------------------------------------------------------
81 // 69 //
82 // Function that returns the intersection of t 70 // Function that returns the intersection of the true path with the surface
83 // of the current volume (either the external 71 // of the current volume (either the external one or the inner one with one
84 // of the daughters: 72 // of the daughters:
85 // 73 //
86 // A = Initial point 74 // A = Initial point
87 // B = another point 75 // B = another point
88 // 76 //
89 // Both A and B are assumed to be on the true 77 // Both A and B are assumed to be on the true path:
90 // 78 //
91 // E is the first point of intersection of 79 // E is the first point of intersection of the chord AB with
92 // a volume other than A (on the surface 80 // a volume other than A (on the surface of A or of a daughter)
93 // 81 //
94 // Convention of Use : 82 // Convention of Use :
95 // i) If it returns "true", then Intersect 83 // i) If it returns "true", then IntersectionPointVelocity is set
96 // to the approximate intersection point 84 // to the approximate intersection point.
97 // ii) If it returns "false", no intersecti 85 // ii) If it returns "false", no intersection was found.
98 // Potential reasons: << 86 // The validity of IntersectedOrRecalculated depends on 'recalculated'
99 // a) no segment found an intersection << 87 // a) if latter is false, then IntersectedOrRecalculated is invalid.
100 // b) too many steps were required - af << 88 // b) if latter is true, then IntersectedOrRecalculated is
101 // and is returning how far it could << 89 // the new endpoint, due to a re-integration.
102 // (If so, it must set 'recalculated <<
103 // TODO/idea: add a new flag: 'unfinish <<
104 // <<
105 // IntersectedOrRecalculated means diff <<
106 // a) if it is the same curve lenght al <<
107 // original enpdoint due to the nee <<
108 // b) if it is at a shorter curve lengt <<
109 // i.e. as far as it could go, becau <<
110 // Note: IntersectedOrRecalculated is v <<
111 // 'true'. <<
112 // ------------------------------------------- 90 // --------------------------------------------------------------------------
113 // NOTE: implementation taken from G4Propagato 91 // NOTE: implementation taken from G4PropagatorInField
114 // 92 //
115 G4bool G4MultiLevelLocator::EstimateIntersecti 93 G4bool G4MultiLevelLocator::EstimateIntersectionPoint(
116 const G4FieldTrack& CurveStart 94 const G4FieldTrack& CurveStartPointVelocity, // A
117 const G4FieldTrack& CurveEndPo 95 const G4FieldTrack& CurveEndPointVelocity, // B
118 const G4ThreeVector& TrialPoint 96 const G4ThreeVector& TrialPoint, // E
119 G4FieldTrack& Intersecte 97 G4FieldTrack& IntersectedOrRecalculatedFT, // Output
120 G4bool& recalculat 98 G4bool& recalculatedEndPoint, // Out
121 G4double& previousSa 99 G4double& previousSafety, // In/Out
122 G4ThreeVector& previousSf 100 G4ThreeVector& previousSftOrigin) // In/Out
123 { 101 {
124 // Find Intersection Point ( A, B, E ) of t 102 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
125 const char* MethodName= "G4MultiLevelLocator << 103
126 <<
127 G4bool found_approximate_intersection = fals 104 G4bool found_approximate_intersection = false;
128 G4bool there_is_no_intersection = fals 105 G4bool there_is_no_intersection = false;
129 << 106
130 G4FieldTrack CurrentA_PointVelocity = Curve 107 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
131 G4FieldTrack CurrentB_PointVelocity = Curve 108 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
132 G4ThreeVector CurrentE_Point = TrialPoint; 109 G4ThreeVector CurrentE_Point = TrialPoint;
133 G4bool validNormalAtE = false; 110 G4bool validNormalAtE = false;
134 G4ThreeVector NormalAtEntry; 111 G4ThreeVector NormalAtEntry;
135 112
136 G4FieldTrack ApproxIntersecPointV(CurveEndP 113 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
137 G4bool validIntersectP= true; // Is <<
138 G4double NewSafety = 0.0; 114 G4double NewSafety = 0.0;
139 G4bool last_AF_intersection = false; 115 G4bool last_AF_intersection = false;
140 116
141 auto integrDriver = GetChordFinderFor()-> << 117 // G4bool final_section= true; // Shows whether current section is last
142 G4bool driverReIntegrates = integrDriver->D << 118 // (i.e. B=full end)
143 <<
144 G4bool first_section = true; 119 G4bool first_section = true;
145 recalculatedEndPoint = false; 120 recalculatedEndPoint = false;
146 G4bool restoredFullEndpoint = false; <<
147 121
148 unsigned int substep_no = 0; << 122 G4bool restoredFullEndpoint = false;
149 123
150 // Statistics for substeps << 124 G4int substep_no = 0;
151 static G4ThreadLocal unsigned int max_no_see <<
152 125
153 #ifdef G4DEBUG_FIELD << 126 G4int oldprc; // cout/cerr precision settings
154 unsigned int trigger_substepno_print = 0; <<
155 const G4double tolerance = 1.0e-8 * CLHEP::m <<
156 unsigned int biggest_depth = 0; <<
157 // using kInitialisingCL = G4LocatorChange <<
158 #endif <<
159 127
160 // Log the location, iteration of changes in << 128 // Limits for substep number
161 //------------------------------------------ << 129 //
162 static G4ThreadLocal G4LocatorChangeLogger e << 130 const G4int max_substeps= 10000; // Test 120 (old value 100 )
163 endChangeB("EndPointB"), recApproxPoint("A << 131 const G4int warn_substeps= 1000; // 100
164 pointH_logger("Trial points 'E': position, <<
165 unsigned int eventCount = 0; <<
166 132
167 if (fCheckMode) << 133 // Statistics for substeps
168 { << 134 //
169 // Clear previous call's data << 135 static G4ThreadLocal G4int max_no_seen= -1;
170 endChangeA.clear(); <<
171 endChangeB.clear(); <<
172 recApproxPoint.clear(); <<
173 pointH_logger.clear(); <<
174 <<
175 // Record the initialisation <<
176 ++eventCount; <<
177 endChangeA.AddRecord( G4LocatorChangeRecor <<
178 eventCount, CurrentA <<
179 endChangeB.AddRecord( G4LocatorChangeRecor <<
180 eventCount, CurrentB <<
181 } <<
182 136
183 //------------------------------------------ 137 //--------------------------------------------------------------------------
184 // Algorithm for the case if progress in fo 138 // Algorithm for the case if progress in founding intersection is too slow.
185 // Process is defined too slow if after N=p 139 // Process is defined too slow if after N=param_substeps advances on the
186 // path, it will be only 'fraction_done' of 140 // path, it will be only 'fraction_done' of the total length.
187 // In this case the remaining length is div 141 // In this case the remaining length is divided in two half and
188 // the loop is restarted for each half. 142 // the loop is restarted for each half.
189 // If progress is still too slow, the divis 143 // If progress is still too slow, the division in two halfs continue
190 // until 'max_depth'. 144 // until 'max_depth'.
191 //------------------------------------------ 145 //--------------------------------------------------------------------------
192 146
193 const G4int param_substeps = 5; // Test val << 147 const G4int param_substeps=5; // Test value for the maximum number
194 // of subst << 148 // of substeps
195 const G4double fraction_done = 0.3; << 149 const G4double fraction_done=0.3;
196 150
197 G4bool Second_half = false; // First half 151 G4bool Second_half = false; // First half or second half of divided step
198 152
199 // We need to know this for the 'final_secti 153 // We need to know this for the 'final_section':
200 // real 'final_section' or first half 'final 154 // real 'final_section' or first half 'final_section'
201 // In algorithm it is considered that the 'S 155 // In algorithm it is considered that the 'Second_half' is true
202 // and it becomes false only if we are in th 156 // and it becomes false only if we are in the first-half of level
203 // depthness or if we are in the first secti 157 // depthness or if we are in the first section
204 158
205 unsigned int depth = 0; // Depth counts subd << 159 G4int depth=0; // Depth counts how many subdivisions of initial step made
206 ++fNumCalls; <<
207 <<
208 NormalAtEntry = GetSurfaceNormal(CurrentE_Po <<
209 160
210 if (fCheckMode) << 161 #ifdef G4DEBUG_FIELD
>> 162 static const G4double tolerance = 1.0e-8 * mm;
>> 163 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
>> 164 if( (TrialPoint - StartPosition).mag() < tolerance)
211 { 165 {
212 pointH_logger.AddRecord( G4LocatorChangeRe << 166 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
213 substep_no, event << 167 "GeomNav1002", JustWarning,
214 G4FieldTrack( Cur << 168 "Intersection point F is exactly at start point A." );
215 0., <<
216 #if (G4DEBUG_FIELD>1) <<
217 G4ThreeVector StartPosition = CurveStartP <<
218 if( (TrialPoint - StartPosition).mag2() < <<
219 { <<
220 ReportImmediateHit( MethodName, StartPo <<
221 tolerance, fNumCall <<
222 } <<
223 #endif <<
224 } 169 }
>> 170 #endif
>> 171
>> 172 NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
225 173
226 // Intermediates Points on the Track = Subdi 174 // Intermediates Points on the Track = Subdivided Points must be stored.
227 // Give the initial values to 'InterMedFt' 175 // Give the initial values to 'InterMedFt'
228 // Important is 'ptrInterMedFT[0]', it saves 176 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
229 // 177 //
230 *ptrInterMedFT[0] = CurveEndPointVelocity; 178 *ptrInterMedFT[0] = CurveEndPointVelocity;
231 for ( auto idepth=1; idepth<max_depth+1; ++i << 179 for (G4int idepth=1; idepth<max_depth+1; idepth++ )
232 { 180 {
233 *ptrInterMedFT[idepth] = CurveStartPointVe << 181 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
234 } 182 }
235 183
236 // Final_section boolean store 184 // Final_section boolean store
237 // 185 //
238 G4bool fin_section_depth[max_depth]; 186 G4bool fin_section_depth[max_depth];
239 for (bool & idepth : fin_section_depth) << 187 for (G4int idepth=0; idepth<max_depth; idepth++ )
240 { 188 {
241 idepth = true; << 189 fin_section_depth[idepth]=true;
242 } 190 }
243 // 'SubStartPoint' is needed to calculate th 191 // 'SubStartPoint' is needed to calculate the length of the divided step
244 // 192 //
245 G4FieldTrack SubStart_PointVelocity = CurveS 193 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
246 194
247 do // Loop checking, 07.10.2016, J.Apostola << 195 do
248 { 196 {
249 unsigned int substep_no_p = 0; << 197 G4int substep_no_p = 0;
250 G4bool sub_final_section = false; // the s 198 G4bool sub_final_section = false; // the same as final_section,
251 // but f 199 // but for 'sub_section'
252 SubStart_PointVelocity = CurrentA_PointVel << 200 SubStart_PointVelocity = CurrentA_PointVelocity;
253 << 201 do // REPEAT param
254 do // Loop checking, 07.10.2016, J.Apostol << 202 {
255 { // REPEAT param <<
256 G4ThreeVector Point_A = CurrentA_PointVe 203 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
257 G4ThreeVector Point_B = CurrentB_PointVe 204 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
258 205
259 #ifdef G4DEBUG_FIELD <<
260 const G4double lenA = CurrentA_PointVelo <<
261 const G4double lenB = CurrentB_PointVelo <<
262 G4double curv_lenAB = lenB - lenA; <<
263 G4double distAB = (Point_B - Point_A <<
264 if( curv_lenAB < distAB * ( 1. - 10.*fiE <<
265 { <<
266 G4cerr << "ERROR> (Start) Point A coin <<
267 << "MLL: iters = " << substep_n <<
268 G4long op=G4cerr.precision(6); <<
269 G4cerr << " Difference = " << di <<
270 << " exceeds limit of relative <<
271 << " i.e. limit = " << 10 * fi <<
272 G4cerr.precision(9); <<
273 G4cerr << " Len A, B = " << len <<
274 << " Position A: " << Po <<
275 << " Position B: " << Po <<
276 G4cerr.precision(op); <<
277 // G4LocatorChangeRecord::ReportVector <<
278 // G4cerr<<"EndPoints A(start) and B(e <<
279 if (fCheckMode) <<
280 { <<
281 G4LocatorChangeLogger::ReportEndChan <<
282 } <<
283 } <<
284 <<
285 if( !validIntersectP ) <<
286 { <<
287 G4ExceptionDescription errmsg; <<
288 errmsg << "Assertion FAILURE - invalid <<
289 << substep_no << " call: " << f <<
290 if (fCheckMode) <<
291 { <<
292 G4LocatorChangeRecord::ReportEndChan <<
293 } <<
294 G4Exception("G4MultiLevelLocator::Esti <<
295 JustWarning, errmsg); <<
296 } <<
297 #endif <<
298 <<
299 // F = a point on true AB path close to 206 // F = a point on true AB path close to point E
300 // (the closest if possible) 207 // (the closest if possible)
301 // 208 //
302 ApproxIntersecPointV = GetChordFinderFor 209 ApproxIntersecPointV = GetChordFinderFor()
303 ->ApproxCurvePoin 210 ->ApproxCurvePointV( CurrentA_PointVelocity,
304 211 CurrentB_PointVelocity,
305 212 CurrentE_Point,
306 << 213 GetEpsilonStepFor());
307 // The above method is the key & most 214 // The above method is the key & most intuitive part ...
308 << 215
309 #ifdef G4DEBUG_FIELD << 216 #ifdef G4DEBUG_FIELD
310 recApproxPoint.push_back(G4LocatorChange << 217 if( ApproxIntersecPointV.GetCurveLength() >
311 << 218 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
312 G4double lenIntsc= ApproxIntersecPointV. <<
313 G4double checkVsEnd= lenB - lenIntsc; <<
314 <<
315 if( lenIntsc > lenB ) <<
316 { 219 {
317 std::ostringstream errmsg; << 220 G4Exception("G4multiLevelLocator::EstimateIntersectionPoint()",
318 errmsg.precision(17); << 221 "GeomNav0003", FatalException,
319 G4double ratio = checkVsEnd / lenB; << 222 "Intermediate F point is past end B point" );
320 G4double ratioTol = std::fabs(ratio) / <<
321 errmsg << "Intermediate F point is pas <<
322 << " l( intersection ) = " << lenInt <<
323 << " l( endpoint ) = " << lenB <<
324 errmsg.precision(8); <<
325 errmsg << " l_end - l_inters = " << <<
326 << " / l_end = " << rati <<
327 << " ratio / tolerance = " << rati <<
328 if( ratioTol < 1.0 ) <<
329 G4Exception(MethodName, "GeomNav0003 <<
330 else <<
331 G4Exception(MethodName, "GeomNav0003 <<
332 } 223 }
333 #endif 224 #endif
334 225
335 G4ThreeVector CurrentF_Point= ApproxInte 226 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
336 227
337 // First check whether EF is small - the 228 // First check whether EF is small - then F is a good approx. point
338 // Calculate the length and direction of 229 // Calculate the length and direction of the chord AF
339 // 230 //
340 G4ThreeVector ChordEF_Vector = CurrentF 231 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
341 232
342 G4ThreeVector NewMomentumDir = ApproxIn << 233 G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
343 G4double MomDir_dot_Norm = NewMome << 234 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
344 << 235
345 #ifdef G4DEBUG_FIELD 236 #ifdef G4DEBUG_FIELD
346 if( fVerboseLevel > 3 ) << 237 if( VerboseLevel > 3 )
347 { 238 {
348 G4ThreeVector ChordAB = Po 239 G4ThreeVector ChordAB = Point_B - Point_A;
349 G4double ABchord_length = Ch 240 G4double ABchord_length = ChordAB.mag();
350 G4double MomDir_dot_ABchord; 241 G4double MomDir_dot_ABchord;
351 MomDir_dot_ABchord = (1.0 / ABchord_l 242 MomDir_dot_ABchord = (1.0 / ABchord_length)
352 * NewMomentumDir.d 243 * NewMomentumDir.dot( ChordAB );
353 G4VIntersectionLocator::ReportTrialSt 244 G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
354 ChordEF_Vector, NewMomentumDir, 245 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
355 G4cout << " dot( MomentumDir, ABchord 246 G4cout << " dot( MomentumDir, ABchord_unit ) = "
356 << MomDir_dot_ABchord << G4end 247 << MomDir_dot_ABchord << G4endl;
357 } 248 }
358 #endif 249 #endif
359 G4bool adequate_angle = 250 G4bool adequate_angle =
360 ( MomDir_dot_Norm >= 0.0 ) // Can 251 ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
361 || (! validNormalAtE ); // Inv 252 || (! validNormalAtE ); // Invalid, cannot use this criterion
362 G4double EF_dist2 = ChordEF_Vector.mag2( 253 G4double EF_dist2 = ChordEF_Vector.mag2();
363 if ( ( EF_dist2 <= sqr(fiDeltaIntersecti 254 if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
364 || ( EF_dist2 <= kCarTolerance*kCarTol 255 || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
365 { 256 {
366 found_approximate_intersection = true; 257 found_approximate_intersection = true;
367 258
368 // Create the "point" return value 259 // Create the "point" return value
369 // 260 //
370 IntersectedOrRecalculatedFT = ApproxIn 261 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
371 IntersectedOrRecalculatedFT.SetPositio 262 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
372 263
373 if ( GetAdjustementOfFoundIntersection 264 if ( GetAdjustementOfFoundIntersection() )
374 { 265 {
375 // Try to Get Correction of Intersec 266 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
376 // 267 //
377 G4ThreeVector IP; 268 G4ThreeVector IP;
378 G4ThreeVector MomentumDir=ApproxInte 269 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
379 G4bool goodCorrection = AdjustmentOf 270 G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
380 CurrentE_P 271 CurrentE_Point, CurrentF_Point, MomentumDir,
381 last_AF_in 272 last_AF_intersection, IP, NewSafety,
382 previousSa 273 previousSafety, previousSftOrigin );
383 if ( goodCorrection ) 274 if ( goodCorrection )
384 { 275 {
385 IntersectedOrRecalculatedFT = Appr 276 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
386 IntersectedOrRecalculatedFT.SetPos 277 IntersectedOrRecalculatedFT.SetPosition(IP);
387 } 278 }
388 } 279 }
389 // Note: in order to return a point on 280 // Note: in order to return a point on the boundary,
390 // we must return E. But it is F 281 // we must return E. But it is F on the curve.
391 // So we must "cheat": we are us 282 // So we must "cheat": we are using the position at point E
392 // and the velocity at point F ! 283 // and the velocity at point F !!!
393 // 284 //
394 // This must limit the length we can a 285 // This must limit the length we can allow for displacement!
395 } 286 }
396 else // E is NOT close enough to the cu 287 else // E is NOT close enough to the curve (ie point F)
397 { 288 {
398 // Check whether any volumes are encou 289 // Check whether any volumes are encountered by the chord AF
399 // ----------------------------------- 290 // ---------------------------------------------------------
400 // First relocate to restore any Voxel 291 // First relocate to restore any Voxel etc information
401 // in the Navigator before calling Com 292 // in the Navigator before calling ComputeStep()
402 // 293 //
403 GetNavigatorFor()->LocateGlobalPointWi 294 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
404 295
405 G4ThreeVector PointG; // Candidate i 296 G4ThreeVector PointG; // Candidate intersection point
406 G4double stepLengthAF; 297 G4double stepLengthAF;
407 G4bool Intersects_FB = false; <<
408 G4bool Intersects_AF = IntersectChord( 298 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
409 299 NewSafety, previousSafety,
410 300 previousSftOrigin,
411 301 stepLengthAF,
412 302 PointG );
413 last_AF_intersection = Intersects_AF; 303 last_AF_intersection = Intersects_AF;
414 if( Intersects_AF ) 304 if( Intersects_AF )
415 { 305 {
416 // G is our new Candidate for the in 306 // G is our new Candidate for the intersection point.
417 // It replaces "E" and we will see << 307 // It replaces "E" and we will repeat the test to see if
418 CurrentB_PointVelocity = ApproxInter << 308 // it is a good enough approximate point for us.
419 CurrentE_Point = PointG; << 309 // B <- F
420 << 310 // E <- G
421 validIntersectP = true; // 'E' ha << 311 //
>> 312 CurrentB_PointVelocity = ApproxIntersecPointV;
>> 313 CurrentE_Point = PointG;
422 314
423 G4bool validNormalLast; 315 G4bool validNormalLast;
424 NormalAtEntry = GetSurfaceNormal( P 316 NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
425 validNormalAtE = validNormalLast; 317 validNormalAtE = validNormalLast;
426 318
427 // As we move point B, must take car << 319 // By moving point B, must take care if current
428 // AF has no intersection to try cur 320 // AF has no intersection to try current FB!!
429 fin_section_depth[depth] = false; << 321 //
430 << 322 fin_section_depth[depth]=false;
431 if (fCheckMode) << 323
432 { << 324
433 ++eventCount; <<
434 endChangeB.push_back( <<
435 G4LocatorChangeRecord(G4LocatorC <<
436 substep_no, event <<
437 } <<
438 #ifdef G4VERBOSE 325 #ifdef G4VERBOSE
439 if( fVerboseLevel > 3 ) 326 if( fVerboseLevel > 3 )
440 { 327 {
441 G4cout << "G4PiF::LI> Investigatin 328 G4cout << "G4PiF::LI> Investigating intermediate point"
442 << " at s=" << ApproxInters 329 << " at s=" << ApproxIntersecPointV.GetCurveLength()
443 << " on way to full s=" 330 << " on way to full s="
444 << CurveEndPointVelocity.Ge 331 << CurveEndPointVelocity.GetCurveLength() << G4endl;
445 } 332 }
446 #endif 333 #endif
447 } 334 }
448 else // not Intersects_AF 335 else // not Intersects_AF
449 { 336 {
450 // In this case: 337 // In this case:
451 // There is NO intersection of AF wi 338 // There is NO intersection of AF with a volume boundary.
452 // We must continue the search in th 339 // We must continue the search in the segment FB!
453 // 340 //
454 GetNavigatorFor()->LocateGlobalPoint 341 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
455 342
456 G4double stepLengthFB; 343 G4double stepLengthFB;
457 G4ThreeVector PointH; 344 G4ThreeVector PointH;
458 345
459 // Check whether any volumes are enc 346 // Check whether any volumes are encountered by the chord FB
460 // --------------------------------- 347 // ---------------------------------------------------------
461 348
462 Intersects_FB = IntersectChord( Curr << 349 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
463 NewS << 350 NewSafety, previousSafety,
464 prev << 351 previousSftOrigin,
465 step << 352 stepLengthFB,
466 Poin << 353 PointH );
467 if( Intersects_FB ) 354 if( Intersects_FB )
468 { 355 {
469 // There is an intersection of FB 356 // There is an intersection of FB with a volume boundary
470 // H <- First Intersection of Chor 357 // H <- First Intersection of Chord FB
471 358
472 // H is our new Candidate for the 359 // H is our new Candidate for the intersection point.
473 // It replaces "E" and we will re 360 // It replaces "E" and we will repeat the test to see if
474 // it is a good enough approximate 361 // it is a good enough approximate point for us.
475 362
476 // Note that F must be in volume v 363 // Note that F must be in volume volA (the same as A)
477 // (otherwise AF would meet a volu 364 // (otherwise AF would meet a volume boundary!)
478 // A <- F 365 // A <- F
479 // E <- H 366 // E <- H
480 // 367 //
481 CurrentA_PointVelocity = ApproxInt 368 CurrentA_PointVelocity = ApproxIntersecPointV;
482 CurrentE_Point = PointH; 369 CurrentE_Point = PointH;
483 370
484 validIntersectP = true; // 'E' <<
485 <<
486 G4bool validNormalH; 371 G4bool validNormalH;
487 NormalAtEntry = GetSurfaceNormal( 372 NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
488 validNormalAtE = validNormalH; << 373 validNormalAtE = validNormalH;
489 374
490 if (fCheckMode) <<
491 { <<
492 ++eventCount; <<
493 endChangeA.push_back( <<
494 G4LocatorChangeRecord(G4Locat <<
495 substep_no, event <<
496 G4FieldTrack intersectH_pn('0'); <<
497 <<
498 intersectH_pn.SetPosition( Point <<
499 intersectH_pn.SetMomentum( Norma <<
500 pointH_logger.AddRecord(G4Locato <<
501 substep_no <<
502 } <<
503 } 375 }
504 else // not Intersects_FB 376 else // not Intersects_FB
505 { 377 {
506 validIntersectP = false; // Int << 378 // There is NO intersection of FB with a volume boundary
507 if( fin_section_depth[depth] ) << 379
>> 380 if(fin_section_depth[depth])
508 { 381 {
509 // If B is the original endpoint 382 // If B is the original endpoint, this means that whatever
510 // volume(s) intersected the ori 383 // volume(s) intersected the original chord, none touch the
511 // smaller chords we have used. 384 // smaller chords we have used.
512 // The value of 'IntersectedOrRe 385 // The value of 'IntersectedOrRecalculatedFT' returned is
513 // likely not valid 386 // likely not valid
514 387
515 // Check on real final_section o 388 // Check on real final_section or SubEndSection
516 // 389 //
517 if( ((Second_half)&&(depth==0)) 390 if( ((Second_half)&&(depth==0)) || (first_section) )
518 { 391 {
519 there_is_no_intersection = tru 392 there_is_no_intersection = true; // real final_section
520 } 393 }
521 else 394 else
522 { 395 {
523 // end of subsection, not real 396 // end of subsection, not real final section
524 // exit from the and go to the 397 // exit from the and go to the depth-1 level
>> 398
525 substep_no_p = param_substeps+ 399 substep_no_p = param_substeps+2; // exit from the loop
526 400
527 // but 'Second_half' is still 401 // but 'Second_half' is still true because we need to find
528 // the 'CurrentE_point' for th 402 // the 'CurrentE_point' for the next loop
>> 403 //
529 Second_half = true; 404 Second_half = true;
530 sub_final_section = true; 405 sub_final_section = true;
>> 406
531 } 407 }
532 } 408 }
533 else 409 else
534 { 410 {
535 CurrentA_PointVelocity = Current << 411 if(depth==0)
536 CurrentB_PointVelocity = (depth= <<
537 <<
538 SubStart_PointVelocity = Current <<
539 restoredFullEndpoint = true; <<
540 <<
541 validIntersectP = false; // 'E' <<
542 <<
543 if (fCheckMode) <<
544 { 412 {
545 ++eventCount; << 413 // We must restore the original endpoint
546 endChangeA.push_back( << 414 //
547 G4LocatorChangeRecord( << 415 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
548 G4LocatorChangeRecord::kNo << 416 CurrentB_PointVelocity = CurveEndPointVelocity;
549 substep_no, event << 417 SubStart_PointVelocity = CurrentA_PointVelocity;
550 endChangeB.push_back( << 418 restoredFullEndpoint = true;
551 G4LocatorChangeRecord ( << 419 }
552 G4LocatorChangeRecord::kNo << 420 else
553 substep_no, event << 421 {
>> 422 // We must restore the depth endpoint
>> 423 //
>> 424 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
>> 425 CurrentB_PointVelocity = *ptrInterMedFT[depth];
>> 426 SubStart_PointVelocity = CurrentA_PointVelocity;
>> 427 restoredFullEndpoint = true;
554 } 428 }
555 } 429 }
556 } // Endif (Intersects_FB) 430 } // Endif (Intersects_FB)
557 } // Endif (Intersects_AF) 431 } // Endif (Intersects_AF)
558 432
559 G4int errorEndPt = 0; // Default: no e << 433 // Ensure that the new endpoints are not further apart in space
>> 434 // than on the curve due to different errors in the integration
>> 435 //
>> 436 G4double linDistSq, curveDist;
>> 437 linDistSq = ( CurrentB_PointVelocity.GetPosition()
>> 438 - CurrentA_PointVelocity.GetPosition() ).mag2();
>> 439 curveDist = CurrentB_PointVelocity.GetCurveLength()
>> 440 - CurrentA_PointVelocity.GetCurveLength();
560 441
561 G4bool recalculatedB= false; << 442 // Change this condition for very strict parameters of propagation
562 if( driverReIntegrates ) << 443 //
563 { << 444 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
564 G4FieldTrack RevisedB_FT = CurrentB_ << 445 {
565 recalculatedB= CheckAndReEstimateEnd << 446 // Re-integrate to obtain a new B
566 << 447 //
567 << 448 G4FieldTrack newEndPointFT=
568 << 449 ReEstimateEndpoint( CurrentA_PointVelocity,
569 if( recalculatedB ) << 450 CurrentB_PointVelocity,
>> 451 linDistSq, // to avoid recalculation
>> 452 curveDist );
>> 453 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
>> 454 CurrentB_PointVelocity = newEndPointFT;
>> 455
>> 456 if ( (fin_section_depth[depth]) // real final section
>> 457 &&( first_section || ((Second_half)&&(depth==0)) ) )
570 { 458 {
571 CurrentB_PointVelocity = RevisedB_ << 459 recalculatedEndPoint = true;
572 // Do not invalidate intersection << 460 IntersectedOrRecalculatedFT = newEndPointFT;
573 // <<
574 // The best course would be to inv <<
575 // BUT if we invalidate it, we mus <<
576 // validIntersectP = false; // <<
577 <<
578 if ( (fin_section_depth[depth]) <<
579 &&( first_section || ((Second_h <<
580 { <<
581 recalculatedEndPoint = true; <<
582 IntersectedOrRecalculatedFT = Re <<
583 // So that we can return it, if 461 // So that we can return it, if it is the endpoint!
584 } <<
585 // else <<
586 // Move forward the other points <<
587 // - or better flag it, so that <<
588 // [ Implementation: a counter <<
589 // => avoids extra work] <<
590 } <<
591 if (fCheckMode) <<
592 { <<
593 ++eventCount; <<
594 endChangeB.push_back( <<
595 G4LocatorChangeRecord( G4Locator <<
596 substep_n <<
597 } 462 }
598 } 463 }
599 else << 464 if( curveDist < 0.0 )
600 { 465 {
601 if( CurrentB_PointVelocity.GetCurveL << 466 fVerboseLevel = 5; // Print out a maximum of information
602 < CurrentA_PointVelocity.GetCurveL << 467 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
603 { << 468 -1.0, NewSafety, substep_no );
604 errorEndPt = 2; << 469 std::ostringstream message;
605 } << 470 message << "Error in advancing propagation." << G4endl
606 } << 471 << " Point A (start) is " << CurrentA_PointVelocity
607 << 472 << G4endl
608 if( errorEndPt > 1 ) // errorEndPt = << 473 << " Point B (end) is " << CurrentB_PointVelocity
609 { << 474 << G4endl
610 std::ostringstream errmsg; << 475 << " Curve distance is " << curveDist << G4endl
>> 476 << G4endl
>> 477 << "The final curve point is not further along"
>> 478 << " than the original!" << G4endl;
611 479
612 ReportReversedPoints(errmsg, << 480 if( recalculatedEndPoint )
613 CurveStartPoint <<
614 NewSafety, fiEp <<
615 CurrentA_PointV <<
616 SubStart_PointV <<
617 ApproxIntersecP <<
618 <<
619 if (fCheckMode) <<
620 { 481 {
621 G4LocatorChangeRecord::ReportEndCh << 482 message << "Recalculation of EndPoint was called with fEpsStep= "
>> 483 << GetEpsilonStepFor() << G4endl;
622 } 484 }
>> 485 oldprc = G4cerr.precision(20);
>> 486 message << " Point A (Curve start) is " << CurveStartPointVelocity
>> 487 << G4endl
>> 488 << " Point B (Curve end) is " << CurveEndPointVelocity
>> 489 << G4endl
>> 490 << " Point A (Current start) is " << CurrentA_PointVelocity
>> 491 << G4endl
>> 492 << " Point B (Current end) is " << CurrentB_PointVelocity
>> 493 << G4endl
>> 494 << " Point S (Sub start) is " << SubStart_PointVelocity
>> 495 << G4endl
>> 496 << " Point E (Trial Point) is " << CurrentE_Point
>> 497 << G4endl
>> 498 << " Point F (Intersection) is " << ApproxIntersecPointV
>> 499 << G4endl
>> 500 << " LocateIntersection parameters are : Substep no= "
>> 501 << substep_no << G4endl
>> 502 << " Substep depth no= "<< substep_no_p << " Depth= "
>> 503 << depth;
>> 504 G4cerr.precision(oldprc);
623 505
624 errmsg << G4endl << " * Location: " << 506 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
625 << "- After EndIf(Intersects_ << 507 "GeomNav0003", FatalException, message);
626 errmsg << " * Bool flags: Recalcula <<
627 << " Intersects_AF = " << I <<
628 << " Intersects_FB = " << I <<
629 errmsg << " * Number of calls to MLL <<
630 G4Exception(MethodName, "GeomNav0003 <<
631 } 508 }
632 if( restoredFullEndpoint ) << 509 if(restoredFullEndpoint)
633 { 510 {
634 fin_section_depth[depth] = restoredF 511 fin_section_depth[depth] = restoredFullEndpoint;
635 restoredFullEndpoint = false; 512 restoredFullEndpoint = false;
636 } 513 }
637 } // EndIf ( E is close enough to the cu 514 } // EndIf ( E is close enough to the curve, ie point F. )
638 // tests ChordAF_Vector.mag() <= maxim 515 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
639 516
640 #ifdef G4DEBUG_FIELD << 517 #ifdef G4DEBUG_LOCATE_INTERSECTION
641 if( trigger_substepno_print == 0) << 518 static G4int trigger_substepno_print= warn_substeps - 20 ;
642 { <<
643 trigger_substepno_print= fWarnSteps - <<
644 } <<
645 519
646 if( substep_no >= trigger_substepno_prin 520 if( substep_no >= trigger_substepno_print )
647 { 521 {
648 G4cout << "Difficulty in converging in << 522 G4cout << "Difficulty in converging in "
>> 523 << "G4MultiLevelLocator::EstimateIntersectionPoint():"
>> 524 << G4endl
649 << " Substep no = " << subst 525 << " Substep no = " << substep_no << G4endl;
650 if( substep_no == trigger_substepno_pr 526 if( substep_no == trigger_substepno_print )
651 { 527 {
652 G4cout << " Start: "; <<
653 printStatus( CurveStartPointVelocity 528 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
654 -1.0, NewSafety, 0 ); << 529 -1.0, NewSafety, 0);
655 if( fCheckMode ) { <<
656 G4LocatorChangeRecord::ReportEndCh <<
657 } else { <<
658 G4cout << " ** For more informatio <<
659 << "-- (it saves and can ou <<
660 } <<
661 } 530 }
662 G4cout << " Point A: "; << 531 G4cout << " State of point A: ";
663 printStatus( CurrentA_PointVelocity, C 532 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
664 -1.0, NewSafety, substep_ << 533 -1.0, NewSafety, substep_no-1, 0);
665 G4cout << " Point B: "; << 534 G4cout << " State of point B: ";
666 printStatus( CurrentA_PointVelocity, C 535 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
667 -1.0, NewSafety, substep_ << 536 -1.0, NewSafety, substep_no);
668 } 537 }
669 #endif 538 #endif
670 ++substep_no; << 539 substep_no++;
671 ++substep_no_p; << 540 substep_no_p++;
672 541
673 } while ( ( ! found_approximate_intersect 542 } while ( ( ! found_approximate_intersection )
674 && ( ! there_is_no_intersection ) 543 && ( ! there_is_no_intersection )
675 && validIntersectP // New c <<
676 && ( substep_no_p <= param_substep 544 && ( substep_no_p <= param_substeps) ); // UNTIL found or
677 545 // failed param substep
>> 546 first_section = false;
678 547
679 if( (!found_approximate_intersection) && ( 548 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
680 { 549 {
681 G4double did_len = std::abs( CurrentA_Po 550 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
682 - SubStart_PointVelocit 551 - SubStart_PointVelocity.GetCurveLength());
683 G4double all_len = std::abs( CurrentB_Po 552 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
684 - SubStart_PointVelocit 553 - SubStart_PointVelocity.GetCurveLength());
685 554
686 G4double distAB = -1; << 555 G4double stepLengthAB;
>> 556 G4ThreeVector PointGe;
>> 557 // Check if progress is too slow and if it possible to go deeper,
>> 558 // then halve the step if so
687 // 559 //
688 // Is progress is too slow, and is it po << 560 if( ( ( did_len )<fraction_done*all_len)
689 // If so, then *halve the step* <<
690 // ============== <<
691 if( (did_len < fraction_done*all_len) <<
692 && (depth<max_depth) && (!sub_final_s 561 && (depth<max_depth) && (!sub_final_section) )
693 { 562 {
694 #ifdef G4DEBUG_FIELD << 563 Second_half=false;
695 static G4ThreadLocal unsigned int numS << 564 depth++;
696 biggest_depth = std::max(depth, bigges <<
697 ++numSplits; <<
698 #endif <<
699 Second_half = false; <<
700 ++depth; <<
701 first_section = false; <<
702 565
703 G4double Sub_len = (all_len-did_len)/( 566 G4double Sub_len = (all_len-did_len)/(2.);
704 G4FieldTrack midPoint = CurrentA_Point << 567 G4FieldTrack start = CurrentA_PointVelocity;
705 G4bool fullAdvance= << 568 G4MagInt_Driver* integrDriver
706 integrDriver->AccurateAdvance(midPo << 569 = GetChordFinderFor()->GetIntegrationDriver();
707 << 570 integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
708 ++fNumAdvanceTrials; << 571 *ptrInterMedFT[depth] = start;
709 if( fullAdvance ) { ++fNumAdvanceFull << 572 CurrentB_PointVelocity = *ptrInterMedFT[depth];
710 << 573
711 G4double lenAchieved= <<
712 midPoint.GetCurveLength()-CurrentA_ <<
713 <<
714 const G4double adequateFraction = (1.0 <<
715 G4bool goodAdvance = (lenAchieved >= a <<
716 if ( goodAdvance ) { ++fNumAdvanceGoo <<
717 <<
718 #ifdef G4DEBUG_FIELD <<
719 else // !goodAdvance <<
720 { <<
721 G4cout << "MLL> AdvanceChordLimited <<
722 << " total length achieved <<
723 << Sub_len << " fraction= " <<
724 if (Sub_len != 0.0 ) { G4cout << le <<
725 else { G4cout << "D <<
726 G4cout << " Good-enough fraction = <<
727 G4cout << " Number of call to mll <<
728 << " iteration " << subste <<
729 << " inner = " << substep_ <<
730 G4cout << " Epsilon of integration <<
731 G4cout << " State at start is <<
732 << G4endl <<
733 << " at end (midpoint <<
734 G4cout << " Particle mass = " << m <<
735 <<
736 G4EquationOfMotion *equation = inte <<
737 ReportFieldValue( CurrentA_PointVel <<
738 ReportFieldValue( midPoint, "midPoi <<
739 G4cout << " Original Start = " <<
740 << CurveStartPointVelocity < <<
741 G4cout << " Original End = " <<
742 << CurveEndPointVelocity << <<
743 G4cout << " Original TrialPoint = <<
744 << TrialPoint << G4endl; <<
745 G4cout << " (this is STRICT mode) <<
746 << " num Splits= " << numSp <<
747 G4cout << G4endl; <<
748 } <<
749 #endif <<
750 <<
751 *ptrInterMedFT[depth] = midPoint; <<
752 CurrentB_PointVelocity = midPoint; <<
753 <<
754 if (fCheckMode) <<
755 { <<
756 ++eventCount; <<
757 endChangeB.push_back( <<
758 G4LocatorChangeRecord( G4LocatorCh <<
759 substep_no, <<
760 } <<
761 <<
762 // Adjust 'SubStartPoint' to calculate 574 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
763 // 575 //
764 SubStart_PointVelocity = CurrentA_Poin 576 SubStart_PointVelocity = CurrentA_PointVelocity;
765 577
766 // Find new trial intersection point n 578 // Find new trial intersection point needed at start of the loop
767 // 579 //
768 G4ThreeVector Point_A = CurrentA_Point 580 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
769 G4ThreeVector Point_B = CurrentB_Point << 581 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
770 << 582
771 G4ThreeVector PointGe; <<
772 GetNavigatorFor()->LocateGlobalPointWi 583 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
773 G4bool Intersects_AB = IntersectChord( << 584 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
774 585 NewSafety, previousSafety,
775 << 586 previousSftOrigin, stepLengthAB,
776 587 PointGe);
777 if( Intersects_AB ) 588 if( Intersects_AB )
778 { 589 {
779 last_AF_intersection = Intersects_AB 590 last_AF_intersection = Intersects_AB;
780 CurrentE_Point = PointGe; 591 CurrentE_Point = PointGe;
781 fin_section_depth[depth] = true; << 592 fin_section_depth[depth]=true;
782 <<
783 validIntersectP = true; // 'E' has <<
784 593
785 G4bool validNormalAB; 594 G4bool validNormalAB;
786 NormalAtEntry = GetSurfaceNormal( P 595 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
787 validNormalAtE = validNormalAB; 596 validNormalAtE = validNormalAB;
788 } 597 }
789 else 598 else
790 { 599 {
791 // No intersection found for first p 600 // No intersection found for first part of curve
792 // (CurrentA,InterMedPoint[depth]). 601 // (CurrentA,InterMedPoint[depth]). Go to the second part
793 // 602 //
794 Second_half = true; 603 Second_half = true;
795 <<
796 validIntersectP= false; // No new <<
797 } 604 }
798 } // if did_len 605 } // if did_len
799 606
800 G4bool unfinished = Second_half; << 607 if( (Second_half)&&(depth!=0) )
801 while ( unfinished && (depth>0) ) // Lo <<
802 { 608 {
803 // Second part of curve (InterMed[dept 609 // Second part of curve (InterMed[depth],Intermed[depth-1]))
804 // On the depth-1 level normally we ar 610 // On the depth-1 level normally we are on the 'second_half'
805 611
>> 612 Second_half = true;
806 // Find new trial intersection point 613 // Find new trial intersection point needed at start of the loop
807 // 614 //
808 SubStart_PointVelocity = *ptrInterMedF 615 SubStart_PointVelocity = *ptrInterMedFT[depth];
809 CurrentA_PointVelocity = *ptrInterMedF 616 CurrentA_PointVelocity = *ptrInterMedFT[depth];
810 CurrentB_PointVelocity = *ptrInterMedF 617 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
811 << 618 // Ensure that the new endpoints are not further apart in space
812 if (fCheckMode) <<
813 { <<
814 ++eventCount; <<
815 G4LocatorChangeRecord chngPop_a( G4L <<
816 substep_no, even <<
817 endChangeA.push_back( chngPop_a ); <<
818 G4LocatorChangeRecord chngPop_b( G4L <<
819 substep_no, even <<
820 endChangeB.push_back( chngPop_b ); <<
821 } <<
822 <<
823 // Ensure that the new endpoints are n <<
824 // than on the curve due to different 619 // than on the curve due to different errors in the integration
825 // 620 //
826 G4int errorEndPt = -1; << 621 G4double linDistSq, curveDist;
827 G4bool recalculatedB= false; << 622 linDistSq = ( CurrentB_PointVelocity.GetPosition()
828 if( driverReIntegrates ) << 623 - CurrentA_PointVelocity.GetPosition() ).mag2();
>> 624 curveDist = CurrentB_PointVelocity.GetCurveLength()
>> 625 - CurrentA_PointVelocity.GetCurveLength();
>> 626 if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
829 { 627 {
830 // Ensure that the new endpoints are << 628 // Re-integrate to obtain a new B
831 // than on the curve due to differen <<
832 // 629 //
833 G4FieldTrack RevisedEndPointFT = Cur << 630 G4FieldTrack newEndPointFT=
834 recalculatedB = << 631 ReEstimateEndpoint( CurrentA_PointVelocity,
835 CheckAndReEstimateEndpoint( Curre << 632 CurrentB_PointVelocity,
836 Curre << 633 linDistSq, // to avoid recalculation
837 Revis << 634 curveDist );
838 error << 635 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
839 if( recalculatedB ) << 636 CurrentB_PointVelocity = newEndPointFT;
840 { << 637 if (depth==1)
841 CurrentB_PointVelocity = RevisedEn << 638 {
842 << 639 recalculatedEndPoint = true;
843 if ( depth == 1 ) << 640 IntersectedOrRecalculatedFT = newEndPointFT;
844 { << 641 // So that we can return it, if it is the endpoint!
845 recalculatedEndPoint = true; <<
846 IntersectedOrRecalculatedFT = R <<
847 // So that we can return it, if <<
848 } <<
849 } 642 }
850 else <<
851 { <<
852 if( CurrentB_PointVelocity.GetCurv <<
853 < CurrentA_PointVelocity.GetCurv <<
854 { <<
855 errorEndPt = 2; <<
856 } <<
857 } <<
858 <<
859 if (fCheckMode) <<
860 { <<
861 ++eventCount; <<
862 endChangeB.push_back( <<
863 G4LocatorChangeRecord(G4LocatorC <<
864 substep_no <<
865 } <<
866 } <<
867 if( errorEndPt > 1 ) // errorEndPt = <<
868 { <<
869 std::ostringstream errmsg; <<
870 ReportReversedPoints(errmsg, <<
871 CurveStartPointVelocity, C <<
872 NewSafety, fiEpsilonStep, <<
873 CurrentA_PointVelocity, Cu <<
874 SubStart_PointVelocity, Cu <<
875 ApproxIntersecPointV, subs <<
876 errmsg << " * Location: " << Method <<
877 errmsg << " * Recalculated = " << re <<
878 G4Exception(MethodName, "GeomNav0003 <<
879 } 643 }
880 << 644
881 G4ThreeVector Point_A = CurrentA_Point << 645 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
882 G4ThreeVector Point_B = CurrentB_Point << 646 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
883 G4ThreeVector PointGi; <<
884 GetNavigatorFor()->LocateGlobalPointWi 647 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
885 G4bool Intersects_AB = IntersectChord( << 648 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
886 649 previousSafety,
887 << 650 previousSftOrigin, stepLengthAB,
888 << 651 PointGe);
889 if( Intersects_AB ) 652 if( Intersects_AB )
890 { 653 {
891 last_AF_intersection = Intersects_AB 654 last_AF_intersection = Intersects_AB;
892 CurrentE_Point = PointGi; << 655 CurrentE_Point = PointGe;
893 656
894 validIntersectP = true; // 'E' has << 657 G4bool validNormalAB;
895 NormalAtEntry = GetSurfaceNormal( P << 658 NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
896 } << 659 validNormalAtE = validNormalAB;
897 else << 660 }
898 { <<
899 validIntersectP = false; // No new <<
900 if( depth == 1) <<
901 { <<
902 there_is_no_intersection = true; <<
903 } <<
904 } <<
905 depth--; 661 depth--;
906 fin_section_depth[depth] = true; << 662 fin_section_depth[depth]=true;
907 unfinished = !validIntersectP; <<
908 } 663 }
909 #ifdef G4DEBUG_FIELD <<
910 if( ! ( validIntersectP || there_is_no_i <<
911 { <<
912 // What happened ?? <<
913 G4cout << "MLL - WARNING Potential FA <<
914 << G4endl <<
915 << " Depth = " << depth << G4e <<
916 << " Num Substeps= " << subste <<
917 G4cout << " Found intersection= " << <<
918 << G4endl; <<
919 G4cout << " Progress report: -- " << <<
920 ReportProgress(G4cout, <<
921 CurveStartPointVelocit <<
922 substep_no, CurrentA_P <<
923 CurrentB_PointVelocity <<
924 NewSafety, depth ); <<
925 G4cout << G4endl; <<
926 } <<
927 #endif <<
928 } // if(!found_aproximate_intersection) 664 } // if(!found_aproximate_intersection)
929 665
930 assert( validIntersectP || there_is_no_int <<
931 || found_approxima <<
932 <<
933 } while ( ( ! found_approximate_intersection 666 } while ( ( ! found_approximate_intersection )
934 && ( ! there_is_no_intersection ) 667 && ( ! there_is_no_intersection )
935 && ( substep_no <= fMaxSteps) ); / << 668 && ( substep_no <= max_substeps) ); // UNTIL found or failed
936 669
937 if( substep_no > max_no_seen ) 670 if( substep_no > max_no_seen )
938 { 671 {
939 max_no_seen = substep_no; 672 max_no_seen = substep_no;
940 #ifdef G4DEBUG_FIELD << 673 #ifdef G4DEBUG_LOCATE_INTERSECTION
941 if( max_no_seen > fWarnSteps ) << 674 if( max_no_seen > warn_substeps )
942 { 675 {
943 trigger_substepno_print = max_no_seen-20 676 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
944 } 677 }
945 #endif 678 #endif
946 } 679 }
947 680
948 if( !there_is_no_intersection && !found_appr << 681 if( ( substep_no >= max_substeps)
>> 682 && !there_is_no_intersection
>> 683 && !found_approximate_intersection )
949 { 684 {
950 if( substep_no >= fMaxSteps) << 685 G4cout << "ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
951 { << 686 << G4endl
952 // Since we cannot go further (yet), w << 687 << " Start and end-point of requested Step:" << G4endl;
953 << 688 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
954 recalculatedEndPoint = true; << 689 -1.0, NewSafety, 0);
955 IntersectedOrRecalculatedFT = CurrentA << 690 G4cout << " Start and end-point of current Sub-Step:" << G4endl;
956 found_approximate_intersection = false << 691 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
957 << 692 -1.0, NewSafety, substep_no-1);
958 std::ostringstream message; << 693 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
959 message << G4endl; << 694 -1.0, NewSafety, substep_no);
960 message << "Convergence is requiring t << 695 std::ostringstream message;
961 << substep_no << " ( limit = << 696 message << "Too many substeps!" << G4endl
962 << G4endl << 697 << " Convergence is requiring too many substeps: "
963 << " Abandoning effort to int << 698 << substep_no << G4endl
964 message << " Number of calls to MLL << 699 << " Abandoning effort to intersect. " << G4endl
965 message << " iteration = " << substep << 700 << " Found intersection = "
966 << 701 << found_approximate_intersection << G4endl
967 message.precision( 12 ); << 702 << " Intersection exists = "
968 G4double done_len = CurrentA_PointVelo << 703 << !there_is_no_intersection << G4endl;
969 G4double full_len = CurveEndPointVeloc << 704
970 message << " Undertaken only le << 705 #ifdef FUTURE_CORRECTION
971 << " out of " << full_len << " << 706 // Attempt to correct the results of the method // FIX - TODO
972 << " Remaining length = <<
973 <<
974 message << " Start and end-point o <<
975 printStatus( CurveStartPointVelocity, <<
976 -1.0, NewSafety, 0, m <<
977 message << " Start and end-poin <<
978 printStatus( CurrentA_PointVelocity, C <<
979 -1.0, NewSafety, substep_ <<
980 printStatus( CurrentA_PointVelocity, C <<
981 -1.0, NewSafety, substep_ <<
982 <<
983 G4Exception(MethodName, "GeomNav0003", <<
984 } <<
985 else if( substep_no >= fWarnSteps) <<
986 { <<
987 std::ostringstream message; <<
988 message << "Many substeps while trying <<
989 << G4endl <<
990 << " Undertaken lengt <<
991 << CurrentB_PointVelocity.GetC <<
992 << " - Needed: " << substep_n <<
993 << " Warning number = <<
994 << " and maximum substeps = " <<
995 G4Exception(MethodName, "GeomNav1002", <<
996 } <<
997 } <<
998 <<
999 return (!there_is_no_intersection) && found <<
1000 // Success or failure <<
1001 } <<
1002 707
1003 void G4MultiLevelLocator::ReportStatistics() << 708 if ( ! found_approximate_intersection )
1004 { << 709 {
1005 G4cout << " Number of calls = " << fNumCal << 710 recalculatedEndPoint = true;
1006 G4cout << " Number of split level ('advanc << 711 // Return the further valid intersection point -- potentially A ??
1007 << fNumAdvanceTrials << G4endl; << 712 // JA/19 Jan 2006
1008 G4cout << " Number of full advances: << 713 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
1009 << fNumAdvanceGood << G4endl; << 714
1010 G4cout << " Number of good advances: << 715 message << " Did not convergence after " << substep_no
1011 << fNumAdvanceFull << G4endl; << 716 << " substeps." << G4endl
1012 } << 717 << " The endpoint was adjused to pointA resulting"
>> 718 << G4endl
>> 719 << " from the last substep: " << CurrentA_PointVelocity
>> 720 << G4endl;
>> 721 }
>> 722 #endif
1013 723
1014 void G4MultiLevelLocator::ReportFieldValue( c << 724 oldprc = G4cout.precision( 10 );
1015 c << 725 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
1016 c << 726 G4double full_len = CurveEndPointVelocity.GetCurveLength();
1017 { << 727 message << " Undertaken only length: " << done_len
1018 enum { maxNumFieldComp = 24 }; << 728 << " out of " << full_len << " required." << G4endl
>> 729 << " Remaining length = " << full_len - done_len;
>> 730 G4cout.precision( oldprc );
1019 731
1020 G4ThreeVector position = locationPV.GetPos << 732 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
1021 G4double startPoint[4] = { position.x(), p << 733 "GeomNav0003", FatalException, message);
1022 locationPV.GetL << 734 }
1023 G4double FieldVec[maxNumFieldComp]; // 24 << 735 else if( substep_no >= warn_substeps )
1024 for (double & i : FieldVec) << 736 {
1025 { << 737 oldprc = G4cout.precision( 10 );
1026 i = 0.0; << 738 std::ostringstream message;
1027 } << 739 message << "Many substeps while trying to locate intersection."
1028 equation->GetFieldValue( startPoint, Field << 740 << G4endl
1029 G4cout << " B-field value (" << nameLoc < << 741 << " Undertaken length: "
1030 << FieldVec[0] << " " << FieldVec[1 << 742 << CurrentB_PointVelocity.GetCurveLength()
1031 G4double Emag2= G4ThreeVector( FieldVec[3] << 743 << " - Needed: " << substep_no << " substeps." << G4endl
1032 FieldVec[4] << 744 << " Warning level = " << warn_substeps
1033 FieldVec[5] << 745 << " and maximum substeps = " << max_substeps;
1034 if( Emag2 > 0.0 ) << 746 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
1035 { << 747 "GeomNav1002", JustWarning, message);
1036 G4cout << " Electric = " << FieldVec[3] << 748 G4cout.precision( oldprc );
1037 << FieldVec[4] << 749 }
1038 << FieldVec[5] << 750 return !there_is_no_intersection; // Success or failure
1039 } <<
1040 return; <<
1041 } 751 }
1042 752