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