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