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