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