| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4MultiLevelLocator.cc,v 1.6 2010-07-13 15:59:42 gcosmo Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
>> 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;
149 122
150 // Statistics for substeps << 123 G4int oldprc; // cout/cerr precision settings
151 static G4ThreadLocal unsigned int max_no_see <<
152 <<
153 #ifdef G4DEBUG_FIELD <<
154 unsigned int trigger_substepno_print = 0; <<
155 const G4double tolerance = 1.0e-8 * CLHEP::m <<
156 unsigned int biggest_depth = 0; <<
157 // using kInitialisingCL = G4LocatorChange <<
158 #endif <<
159 124
160 // Log the location, iteration of changes in << 125 // Limits for substep number
161 //------------------------------------------ << 126 //
162 static G4ThreadLocal G4LocatorChangeLogger e << 127 const G4int max_substeps= 10000; // Test 120 (old value 100 )
163 endChangeB("EndPointB"), recApproxPoint("A << 128 const G4int warn_substeps= 1000; // 100
164 pointH_logger("Trial points 'E': position, <<
165 unsigned int eventCount = 0; <<
166 129
167 if (fCheckMode) << 130 // Statistics for substeps
168 { << 131 //
169 // Clear previous call's data << 132 static G4int max_no_seen= -1;
170 endChangeA.clear(); << 133 static G4int trigger_substepno_print= warn_substeps - 20 ;
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 134
183 //------------------------------------------ 135 //--------------------------------------------------------------------------
184 // Algorithm for the case if progress in fo 136 // Algorithm for the case if progress in founding intersection is too slow.
185 // Process is defined too slow if after N=p 137 // Process is defined too slow if after N=param_substeps advances on the
186 // path, it will be only 'fraction_done' of 138 // path, it will be only 'fraction_done' of the total length.
187 // In this case the remaining length is div 139 // In this case the remaining length is divided in two half and
188 // the loop is restarted for each half. 140 // the loop is restarted for each half.
189 // If progress is still too slow, the divis 141 // If progress is still too slow, the division in two halfs continue
190 // until 'max_depth'. 142 // until 'max_depth'.
191 //------------------------------------------ 143 //--------------------------------------------------------------------------
192 144
193 const G4int param_substeps = 5; // Test val << 145 const G4int param_substeps=10; // Test value for the maximum number
194 // of subst << 146 // of substeps
195 const G4double fraction_done = 0.3; << 147 const G4double fraction_done=0.3;
196 148
197 G4bool Second_half = false; // First half 149 G4bool Second_half = false; // First half or second half of divided step
198 150
199 // We need to know this for the 'final_secti 151 // We need to know this for the 'final_section':
200 // real 'final_section' or first half 'final 152 // real 'final_section' or first half 'final_section'
201 // In algorithm it is considered that the 'S 153 // In algorithm it is considered that the 'Second_half' is true
202 // and it becomes false only if we are in th 154 // and it becomes false only if we are in the first-half of level
203 // depthness or if we are in the first secti 155 // depthness or if we are in the first section
204 156
205 unsigned int depth = 0; // Depth counts subd << 157 G4int depth=0; // Depth counts how many subdivisions of initial step made
206 ++fNumCalls; <<
207 <<
208 NormalAtEntry = GetSurfaceNormal(CurrentE_Po <<
209 158
210 if (fCheckMode) << 159 #ifdef G4DEBUG_FIELD
>> 160 static const G4double tolerance = 1.0e-8 * mm;
>> 161 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
>> 162 if( (TrialPoint - StartPosition).mag() < tolerance)
211 { 163 {
212 pointH_logger.AddRecord( G4LocatorChangeRe << 164 G4cerr << "WARNING - G4MultiLevelLocator::EstimateIntersectionPoint()"
213 substep_no, event << 165 << G4endl
214 G4FieldTrack( Cur << 166 << " Intermediate F point is on top of starting point A."
215 0., << 167 << G4endl;
216 #if (G4DEBUG_FIELD>1) << 168 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
217 G4ThreeVector StartPosition = CurveStartP << 169 "IntersectionPointIsAtStart", JustWarning,
218 if( (TrialPoint - StartPosition).mag2() < << 170 "Intersection point F is exactly at start point A." );
219 { <<
220 ReportImmediateHit( MethodName, StartPo <<
221 tolerance, fNumCall <<
222 } <<
223 #endif <<
224 } 171 }
>> 172 #endif
225 173
226 // Intermediates Points on the Track = Subdi 174 // Intermediates Points on the Track = Subdivided Points must be stored.
227 // Give the initial values to 'InterMedFt' 175 // Give the initial values to 'InterMedFt'
228 // Important is 'ptrInterMedFT[0]', it saves 176 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
229 // 177 //
230 *ptrInterMedFT[0] = CurveEndPointVelocity; 178 *ptrInterMedFT[0] = CurveEndPointVelocity;
231 for ( auto idepth=1; idepth<max_depth+1; ++i << 179 for (G4int idepth=1; idepth<max_depth+1; idepth++ )
232 { 180 {
233 *ptrInterMedFT[idepth] = CurveStartPointVe << 181 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
234 } 182 }
235 183
236 // Final_section boolean store 184 // Final_section boolean store
237 // 185 //
238 G4bool fin_section_depth[max_depth]; 186 G4bool fin_section_depth[max_depth];
239 for (bool & idepth : fin_section_depth) << 187 for (G4int idepth=0; idepth<max_depth; idepth++ )
240 { 188 {
241 idepth = true; << 189 fin_section_depth[idepth]=true;
242 } 190 }
243 // 'SubStartPoint' is needed to calculate th 191 // 'SubStartPoint' is needed to calculate the length of the divided step
244 // 192 //
245 G4FieldTrack SubStart_PointVelocity = CurveS 193 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
246 194
247 do // Loop checking, 07.10.2016, J.Apostola << 195 do
248 { 196 {
249 unsigned int substep_no_p = 0; << 197 G4int substep_no_p = 0;
250 G4bool sub_final_section = false; // the s 198 G4bool sub_final_section = false; // the same as final_section,
251 // but f 199 // but for 'sub_section'
252 SubStart_PointVelocity = CurrentA_PointVel << 200 SubStart_PointVelocity = CurrentA_PointVelocity;
253 << 201 do // REPEAT param
254 do // Loop checking, 07.10.2016, J.Apostol << 202 {
255 { // REPEAT param <<
256 G4ThreeVector Point_A = CurrentA_PointVe 203 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
257 G4ThreeVector Point_B = CurrentB_PointVe 204 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
258 205
259 #ifdef G4DEBUG_FIELD <<
260 const G4double lenA = CurrentA_PointVelo <<
261 const G4double lenB = CurrentB_PointVelo <<
262 G4double curv_lenAB = lenB - lenA; <<
263 G4double distAB = (Point_B - Point_A <<
264 if( curv_lenAB < distAB * ( 1. - 10.*fiE <<
265 { <<
266 G4cerr << "ERROR> (Start) Point A coin <<
267 << "MLL: iters = " << substep_n <<
268 G4long op=G4cerr.precision(6); <<
269 G4cerr << " Difference = " << di <<
270 << " exceeds limit of relative <<
271 << " i.e. limit = " << 10 * fi <<
272 G4cerr.precision(9); <<
273 G4cerr << " Len A, B = " << len <<
274 << " Position A: " << Po <<
275 << " Position B: " << Po <<
276 G4cerr.precision(op); <<
277 // G4LocatorChangeRecord::ReportVector <<
278 // G4cerr<<"EndPoints A(start) and B(e <<
279 if (fCheckMode) <<
280 { <<
281 G4LocatorChangeLogger::ReportEndChan <<
282 } <<
283 } <<
284 <<
285 if( !validIntersectP ) <<
286 { <<
287 G4ExceptionDescription errmsg; <<
288 errmsg << "Assertion FAILURE - invalid <<
289 << substep_no << " call: " << f <<
290 if (fCheckMode) <<
291 { <<
292 G4LocatorChangeRecord::ReportEndChan <<
293 } <<
294 G4Exception("G4MultiLevelLocator::Esti <<
295 JustWarning, errmsg); <<
296 } <<
297 #endif <<
298 <<
299 // F = a point on true AB path close to 206 // F = a point on true AB path close to point E
300 // (the closest if possible) 207 // (the closest if possible)
301 // 208 //
302 ApproxIntersecPointV = GetChordFinderFor 209 ApproxIntersecPointV = GetChordFinderFor()
303 ->ApproxCurvePoin 210 ->ApproxCurvePointV( CurrentA_PointVelocity,
304 211 CurrentB_PointVelocity,
305 212 CurrentE_Point,
306 << 213 GetEpsilonStepFor());
307 // The above method is the key & most 214 // The above method is the key & most intuitive part ...
308 << 215
309 #ifdef G4DEBUG_FIELD << 216 #ifdef G4DEBUG_FIELD
310 recApproxPoint.push_back(G4LocatorChange << 217 if( ApproxIntersecPointV.GetCurveLength() >
311 << 218 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
312 G4double lenIntsc= ApproxIntersecPointV. <<
313 G4double checkVsEnd= lenB - lenIntsc; <<
314 <<
315 if( lenIntsc > lenB ) <<
316 { 219 {
317 std::ostringstream errmsg; << 220 G4cerr << "ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
318 errmsg.precision(17); << 221 << G4endl
319 G4double ratio = checkVsEnd / lenB; << 222 << " Intermediate F point is more advanced than"
320 G4double ratioTol = std::fabs(ratio) / << 223 << " endpoint B." << G4endl;
321 errmsg << "Intermediate F point is pas << 224 G4Exception("G4multiLevelLocator::EstimateIntersectionPoint()",
322 << " l( intersection ) = " << lenInt << 225 "IntermediatePointConfusion", FatalException,
323 << " l( endpoint ) = " << lenB << 226 "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 } 227 }
333 #endif 228 #endif
334 229
335 G4ThreeVector CurrentF_Point= ApproxInte 230 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
336 231
337 // First check whether EF is small - the 232 // First check whether EF is small - then F is a good approx. point
338 // Calculate the length and direction of 233 // Calculate the length and direction of the chord AF
339 // 234 //
340 G4ThreeVector ChordEF_Vector = CurrentF 235 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
341 236
342 G4ThreeVector NewMomentumDir = ApproxIn << 237 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 { 238 {
366 found_approximate_intersection = true; 239 found_approximate_intersection = true;
367 240
368 // Create the "point" return value 241 // Create the "point" return value
369 // 242 //
370 IntersectedOrRecalculatedFT = ApproxIn 243 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
371 IntersectedOrRecalculatedFT.SetPositio 244 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
372 245
373 if ( GetAdjustementOfFoundIntersection 246 if ( GetAdjustementOfFoundIntersection() )
374 { 247 {
375 // Try to Get Correction of Intersec 248 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
376 // 249 //
377 G4ThreeVector IP; 250 G4ThreeVector IP;
378 G4ThreeVector MomentumDir=ApproxInte 251 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
379 G4bool goodCorrection = AdjustmentOf 252 G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
380 CurrentE_P 253 CurrentE_Point, CurrentF_Point, MomentumDir,
381 last_AF_in 254 last_AF_intersection, IP, NewSafety,
382 previousSa << 255 fPreviousSafety, fPreviousSftOrigin );
383 if ( goodCorrection ) 256 if ( goodCorrection )
384 { 257 {
385 IntersectedOrRecalculatedFT = Appr 258 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
386 IntersectedOrRecalculatedFT.SetPos 259 IntersectedOrRecalculatedFT.SetPosition(IP);
387 } 260 }
388 } 261 }
389 // Note: in order to return a point on 262 // Note: in order to return a point on the boundary,
390 // we must return E. But it is F 263 // we must return E. But it is F on the curve.
391 // So we must "cheat": we are us 264 // So we must "cheat": we are using the position at point E
392 // and the velocity at point F ! 265 // and the velocity at point F !!!
393 // 266 //
394 // This must limit the length we can a 267 // This must limit the length we can allow for displacement!
395 } 268 }
396 else // E is NOT close enough to the cu 269 else // E is NOT close enough to the curve (ie point F)
397 { 270 {
398 // Check whether any volumes are encou 271 // Check whether any volumes are encountered by the chord AF
399 // ----------------------------------- 272 // ---------------------------------------------------------
400 // First relocate to restore any Voxel 273 // First relocate to restore any Voxel etc information
401 // in the Navigator before calling Com 274 // in the Navigator before calling ComputeStep()
402 // 275 //
403 GetNavigatorFor()->LocateGlobalPointWi 276 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
404 277
405 G4ThreeVector PointG; // Candidate i 278 G4ThreeVector PointG; // Candidate intersection point
406 G4double stepLengthAF; 279 G4double stepLengthAF;
407 G4bool Intersects_FB = false; <<
408 G4bool Intersects_AF = IntersectChord( 280 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
409 << 281 NewSafety,fPreviousSafety,
410 << 282 fPreviousSftOrigin,
411 283 stepLengthAF,
412 284 PointG );
413 last_AF_intersection = Intersects_AF; 285 last_AF_intersection = Intersects_AF;
414 if( Intersects_AF ) 286 if( Intersects_AF )
415 { 287 {
416 // G is our new Candidate for the in 288 // G is our new Candidate for the intersection point.
417 // It replaces "E" and we will see << 289 // It replaces "E" and we will repeat the test to see if
418 CurrentB_PointVelocity = ApproxInter << 290 // it is a good enough approximate point for us.
419 CurrentE_Point = PointG; << 291 // B <- F
420 << 292 // E <- G
421 validIntersectP = true; // 'E' ha << 293 //
422 << 294 CurrentB_PointVelocity = ApproxIntersecPointV;
423 G4bool validNormalLast; << 295 CurrentE_Point = PointG;
424 NormalAtEntry = GetSurfaceNormal( P <<
425 validNormalAtE = validNormalLast; <<
426 296
427 // As we move point B, must take car << 297 // By moving point B, must take care if current
428 // AF has no intersection to try cur 298 // AF has no intersection to try current FB!!
429 fin_section_depth[depth] = false; << 299 //
430 << 300 fin_section_depth[depth]=false;
431 if (fCheckMode) << 301
432 { << 302
433 ++eventCount; <<
434 endChangeB.push_back( <<
435 G4LocatorChangeRecord(G4LocatorC <<
436 substep_no, event <<
437 } <<
438 #ifdef G4VERBOSE 303 #ifdef G4VERBOSE
439 if( fVerboseLevel > 3 ) 304 if( fVerboseLevel > 3 )
440 { 305 {
441 G4cout << "G4PiF::LI> Investigatin 306 G4cout << "G4PiF::LI> Investigating intermediate point"
442 << " at s=" << ApproxInters 307 << " at s=" << ApproxIntersecPointV.GetCurveLength()
443 << " on way to full s=" 308 << " on way to full s="
444 << CurveEndPointVelocity.Ge 309 << CurveEndPointVelocity.GetCurveLength() << G4endl;
445 } 310 }
446 #endif 311 #endif
447 } 312 }
448 else // not Intersects_AF 313 else // not Intersects_AF
449 { 314 {
450 // In this case: 315 // In this case:
451 // There is NO intersection of AF wi 316 // There is NO intersection of AF with a volume boundary.
452 // We must continue the search in th 317 // We must continue the search in the segment FB!
453 // 318 //
454 GetNavigatorFor()->LocateGlobalPoint 319 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
455 320
456 G4double stepLengthFB; 321 G4double stepLengthFB;
457 G4ThreeVector PointH; 322 G4ThreeVector PointH;
458 323
459 // Check whether any volumes are enc 324 // Check whether any volumes are encountered by the chord FB
460 // --------------------------------- 325 // ---------------------------------------------------------
461 326
462 Intersects_FB = IntersectChord( Curr << 327 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
463 NewS << 328 NewSafety,fPreviousSafety,
464 prev << 329 fPreviousSftOrigin,
465 step << 330 stepLengthFB,
466 Poin << 331 PointH );
467 if( Intersects_FB ) 332 if( Intersects_FB )
468 { 333 {
469 // There is an intersection of FB 334 // There is an intersection of FB with a volume boundary
470 // H <- First Intersection of Chor 335 // H <- First Intersection of Chord FB
471 336
472 // H is our new Candidate for the 337 // H is our new Candidate for the intersection point.
473 // It replaces "E" and we will re 338 // It replaces "E" and we will repeat the test to see if
474 // it is a good enough approximate 339 // it is a good enough approximate point for us.
475 340
476 // Note that F must be in volume v 341 // Note that F must be in volume volA (the same as A)
477 // (otherwise AF would meet a volu 342 // (otherwise AF would meet a volume boundary!)
478 // A <- F 343 // A <- F
479 // E <- H 344 // E <- H
480 // 345 //
481 CurrentA_PointVelocity = ApproxInt 346 CurrentA_PointVelocity = ApproxIntersecPointV;
482 CurrentE_Point = PointH; 347 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 } 348 }
504 else // not Intersects_FB 349 else // not Intersects_FB
505 { 350 {
506 validIntersectP = false; // Int << 351 // There is NO intersection of FB with a volume boundary
507 if( fin_section_depth[depth] ) << 352
>> 353 if(fin_section_depth[depth])
508 { 354 {
509 // If B is the original endpoint 355 // If B is the original endpoint, this means that whatever
510 // volume(s) intersected the ori 356 // volume(s) intersected the original chord, none touch the
511 // smaller chords we have used. 357 // smaller chords we have used.
512 // The value of 'IntersectedOrRe 358 // The value of 'IntersectedOrRecalculatedFT' returned is
513 // likely not valid 359 // likely not valid
514 360
515 // Check on real final_section o 361 // Check on real final_section or SubEndSection
516 // 362 //
517 if( ((Second_half)&&(depth==0)) 363 if( ((Second_half)&&(depth==0)) || (first_section) )
518 { 364 {
519 there_is_no_intersection = tru 365 there_is_no_intersection = true; // real final_section
520 } 366 }
521 else 367 else
522 { 368 {
523 // end of subsection, not real 369 // end of subsection, not real final section
524 // exit from the and go to the 370 // exit from the and go to the depth-1 level
>> 371
525 substep_no_p = param_substeps+ 372 substep_no_p = param_substeps+2; // exit from the loop
526 373
527 // but 'Second_half' is still 374 // but 'Second_half' is still true because we need to find
528 // the 'CurrentE_point' for th 375 // the 'CurrentE_point' for the next loop
>> 376 //
529 Second_half = true; 377 Second_half = true;
530 sub_final_section = true; 378 sub_final_section = true;
>> 379
531 } 380 }
532 } 381 }
533 else 382 else
534 { 383 {
535 CurrentA_PointVelocity = Current << 384 if(depth==0)
536 CurrentB_PointVelocity = (depth= << 385 {
537 << 386 // We must restore the original endpoint
538 SubStart_PointVelocity = Current << 387 //
539 restoredFullEndpoint = true; << 388 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
540 << 389 CurrentB_PointVelocity = CurveEndPointVelocity;
541 validIntersectP = false; // 'E' << 390 SubStart_PointVelocity = CurrentA_PointVelocity;
542 << 391 restoredFullEndpoint = true;
543 if (fCheckMode) << 392 }
>> 393 else
544 { 394 {
545 ++eventCount; << 395 // We must restore the depth endpoint
546 endChangeA.push_back( << 396 //
547 G4LocatorChangeRecord( << 397 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
548 G4LocatorChangeRecord::kNo << 398 CurrentB_PointVelocity = *ptrInterMedFT[depth];
549 substep_no, event << 399 SubStart_PointVelocity = CurrentA_PointVelocity;
550 endChangeB.push_back( << 400 restoredFullEndpoint = true;
551 G4LocatorChangeRecord ( <<
552 G4LocatorChangeRecord::kNo <<
553 substep_no, event <<
554 } 401 }
555 } 402 }
556 } // Endif (Intersects_FB) 403 } // Endif (Intersects_FB)
557 } // Endif (Intersects_AF) 404 } // Endif (Intersects_AF)
558 405
559 G4int errorEndPt = 0; // Default: no e << 406 // Ensure that the new endpoints are not further apart in space
>> 407 // than on the curve due to different errors in the integration
>> 408 //
>> 409 G4double linDistSq, curveDist;
>> 410 linDistSq = ( CurrentB_PointVelocity.GetPosition()
>> 411 - CurrentA_PointVelocity.GetPosition() ).mag2();
>> 412 curveDist = CurrentB_PointVelocity.GetCurveLength()
>> 413 - CurrentA_PointVelocity.GetCurveLength();
560 414
561 G4bool recalculatedB= false; << 415 // Change this condition for very strict parameters of propagation
562 if( driverReIntegrates ) << 416 //
563 { << 417 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
564 G4FieldTrack RevisedB_FT = CurrentB_ << 418 {
565 recalculatedB= CheckAndReEstimateEnd << 419 // Re-integrate to obtain a new B
566 << 420 //
567 << 421 G4FieldTrack newEndPointFT=
568 << 422 ReEstimateEndpoint( CurrentA_PointVelocity,
569 if( recalculatedB ) << 423 CurrentB_PointVelocity,
>> 424 linDistSq, // to avoid recalculation
>> 425 curveDist );
>> 426 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
>> 427 CurrentB_PointVelocity = newEndPointFT;
>> 428
>> 429 if ( (fin_section_depth[depth]) // real final section
>> 430 &&( first_section || ((Second_half)&&(depth==0)) ) )
570 { 431 {
571 CurrentB_PointVelocity = RevisedB_ << 432 recalculatedEndPoint = true;
572 // Do not invalidate intersection << 433 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 434 // 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 } 435 }
598 } 436 }
599 else << 437 if( curveDist < 0.0 )
600 { 438 {
601 if( CurrentB_PointVelocity.GetCurveL << 439 G4cerr << "ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
602 < CurrentA_PointVelocity.GetCurveL << 440 << G4endl
603 { << 441 << " Error in advancing propagation." << G4endl;
604 errorEndPt = 2; << 442 fVerboseLevel = 5; // Print out a maximum of information
605 } << 443 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
>> 444 -1.0, NewSafety, substep_no );
>> 445 G4cerr << " Point A (start) is " << CurrentA_PointVelocity
>> 446 << G4endl;
>> 447 G4cerr << " Point B (end) is " << CurrentB_PointVelocity
>> 448 << G4endl;
>> 449 G4cerr << " Curve distance is " << curveDist << G4endl;
>> 450 G4cerr << G4endl
>> 451 << "The final curve point is not further along"
>> 452 << " than the original!" << G4endl;
>> 453
>> 454 if( recalculatedEndPoint )
>> 455 {
>> 456 G4cerr << "Recalculation of EndPoint was called with fEpsStep= "
>> 457 << GetEpsilonStepFor() << G4endl;
>> 458 }
>> 459 oldprc = G4cerr.precision(20);
>> 460 G4cerr << " Point A (Curve start) is " << CurveStartPointVelocity
>> 461 << G4endl;
>> 462 G4cerr << " Point B (Curve end) is " << CurveEndPointVelocity
>> 463 << G4endl;
>> 464 G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity
>> 465 << G4endl;
>> 466 G4cerr << " Point B (Current end) is " << CurrentB_PointVelocity
>> 467 << G4endl;
>> 468 G4cerr << " Point S (Sub start) is " << SubStart_PointVelocity
>> 469 << G4endl;
>> 470 G4cerr << " Point E (Trial Point) is " << CurrentE_Point
>> 471 << G4endl;
>> 472 G4cerr << " Point F (Intersection) is " << ApproxIntersecPointV
>> 473 << G4endl;
>> 474 G4cerr << " LocateIntersection parameters are : Substep no= "
>> 475 << substep_no << G4endl;
>> 476 G4cerr << " Substep depth no= "<< substep_no_p << " Depth= "
>> 477 << depth << G4endl;
>> 478 G4cerr.precision(oldprc);
>> 479
>> 480 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
>> 481 "FatalError", FatalException,
>> 482 "Error in advancing propagation.");
606 } 483 }
607 << 484 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 { 485 {
634 fin_section_depth[depth] = restoredF 486 fin_section_depth[depth] = restoredFullEndpoint;
635 restoredFullEndpoint = false; 487 restoredFullEndpoint = false;
636 } 488 }
637 } // EndIf ( E is close enough to the cu 489 } // EndIf ( E is close enough to the curve, ie point F. )
638 // tests ChordAF_Vector.mag() <= maxim 490 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
639 491
640 #ifdef G4DEBUG_FIELD << 492 #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 493 if( substep_no >= trigger_substepno_print )
647 { 494 {
648 G4cout << "Difficulty in converging in << 495 G4cout << "Difficulty in converging in "
>> 496 << "G4MultiLevelLocator::EstimateIntersectionPoint():"
>> 497 << G4endl
649 << " Substep no = " << subst 498 << " Substep no = " << substep_no << G4endl;
650 if( substep_no == trigger_substepno_pr 499 if( substep_no == trigger_substepno_print )
651 { 500 {
652 G4cout << " Start: "; <<
653 printStatus( CurveStartPointVelocity 501 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
654 -1.0, NewSafety, 0 ); << 502 -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 } 503 }
662 G4cout << " Point A: "; << 504 G4cout << " State of point A: ";
663 printStatus( CurrentA_PointVelocity, C 505 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
664 -1.0, NewSafety, substep_ << 506 -1.0, NewSafety, substep_no-1, 0);
665 G4cout << " Point B: "; << 507 G4cout << " State of point B: ";
666 printStatus( CurrentA_PointVelocity, C 508 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
667 -1.0, NewSafety, substep_ << 509 -1.0, NewSafety, substep_no);
668 } 510 }
669 #endif 511 #endif
670 ++substep_no; << 512 substep_no++;
671 ++substep_no_p; << 513 substep_no_p++;
672 514
673 } while ( ( ! found_approximate_intersect 515 } while ( ( ! found_approximate_intersection )
674 && ( ! there_is_no_intersection ) 516 && ( ! there_is_no_intersection )
675 && validIntersectP // New c <<
676 && ( substep_no_p <= param_substep 517 && ( substep_no_p <= param_substeps) ); // UNTIL found or
677 518 // failed param substep
>> 519 first_section = false;
678 520
679 if( (!found_approximate_intersection) && ( 521 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
680 { 522 {
681 G4double did_len = std::abs( CurrentA_Po 523 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
682 - SubStart_PointVelocit 524 - SubStart_PointVelocity.GetCurveLength());
683 G4double all_len = std::abs( CurrentB_Po 525 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
684 - SubStart_PointVelocit 526 - SubStart_PointVelocity.GetCurveLength());
685 527
686 G4double distAB = -1; << 528 G4double stepLengthAB;
>> 529 G4ThreeVector PointGe;
>> 530 // Check if progress is too slow and if it possible to go deeper,
>> 531 // then halve the step if so
687 // 532 //
688 // Is progress is too slow, and is it po << 533 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 534 && (depth<max_depth) && (!sub_final_section) )
693 { 535 {
694 #ifdef G4DEBUG_FIELD << 536 Second_half=false;
695 static G4ThreadLocal unsigned int numS << 537 depth++;
696 biggest_depth = std::max(depth, bigges <<
697 ++numSplits; <<
698 #endif <<
699 Second_half = false; <<
700 ++depth; <<
701 first_section = false; <<
702 538
703 G4double Sub_len = (all_len-did_len)/( 539 G4double Sub_len = (all_len-did_len)/(2.);
704 G4FieldTrack midPoint = CurrentA_Point << 540 G4FieldTrack start = CurrentA_PointVelocity;
705 G4bool fullAdvance= << 541 G4MagInt_Driver* integrDriver
706 integrDriver->AccurateAdvance(midPo << 542 = GetChordFinderFor()->GetIntegrationDriver();
707 << 543 integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
708 ++fNumAdvanceTrials; << 544 *ptrInterMedFT[depth] = start;
709 if( fullAdvance ) { ++fNumAdvanceFull << 545 CurrentB_PointVelocity = *ptrInterMedFT[depth];
710 << 546
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 547 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
763 // 548 //
764 SubStart_PointVelocity = CurrentA_Poin 549 SubStart_PointVelocity = CurrentA_PointVelocity;
765 550
766 // Find new trial intersection point n 551 // Find new trial intersection point needed at start of the loop
767 // 552 //
768 G4ThreeVector Point_A = CurrentA_Point 553 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
769 G4ThreeVector Point_B = CurrentB_Point << 554 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
770 << 555
771 G4ThreeVector PointGe; <<
772 GetNavigatorFor()->LocateGlobalPointWi 556 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
773 G4bool Intersects_AB = IntersectChord( << 557 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
774 << 558 NewSafety, fPreviousSafety,
775 << 559 fPreviousSftOrigin,stepLengthAB,
776 560 PointGe);
777 if( Intersects_AB ) 561 if( Intersects_AB )
778 { 562 {
779 last_AF_intersection = Intersects_AB 563 last_AF_intersection = Intersects_AB;
780 CurrentE_Point = PointGe; 564 CurrentE_Point = PointGe;
781 fin_section_depth[depth] = true; << 565 fin_section_depth[depth]=true;
782 <<
783 validIntersectP = true; // 'E' has <<
784 <<
785 G4bool validNormalAB; <<
786 NormalAtEntry = GetSurfaceNormal( P <<
787 validNormalAtE = validNormalAB; <<
788 } 566 }
789 else 567 else
790 { 568 {
791 // No intersection found for first p 569 // No intersection found for first part of curve
792 // (CurrentA,InterMedPoint[depth]). 570 // (CurrentA,InterMedPoint[depth]). Go to the second part
793 // 571 //
794 Second_half = true; 572 Second_half = true;
795 <<
796 validIntersectP= false; // No new <<
797 } 573 }
798 } // if did_len 574 } // if did_len
799 575
800 G4bool unfinished = Second_half; << 576 if( (Second_half)&&(depth!=0) )
801 while ( unfinished && (depth>0) ) // Lo <<
802 { 577 {
803 // Second part of curve (InterMed[dept << 578 // Second part of curve (InterMed[depth],Intermed[depth-1]) )
804 // On the depth-1 level normally we ar 579 // On the depth-1 level normally we are on the 'second_half'
805 580
>> 581 Second_half = true;
806 // Find new trial intersection point 582 // Find new trial intersection point needed at start of the loop
807 // 583 //
808 SubStart_PointVelocity = *ptrInterMedF 584 SubStart_PointVelocity = *ptrInterMedFT[depth];
809 CurrentA_PointVelocity = *ptrInterMedF 585 CurrentA_PointVelocity = *ptrInterMedFT[depth];
810 CurrentB_PointVelocity = *ptrInterMedF 586 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
811 << 587 // 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 588 // than on the curve due to different errors in the integration
825 // 589 //
826 G4int errorEndPt = -1; << 590 G4double linDistSq, curveDist;
827 G4bool recalculatedB= false; << 591 linDistSq = ( CurrentB_PointVelocity.GetPosition()
828 if( driverReIntegrates ) << 592 - CurrentA_PointVelocity.GetPosition() ).mag2();
>> 593 curveDist = CurrentB_PointVelocity.GetCurveLength()
>> 594 - CurrentA_PointVelocity.GetCurveLength();
>> 595 if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
829 { 596 {
830 // Ensure that the new endpoints are << 597 // Re-integrate to obtain a new B
831 // than on the curve due to differen <<
832 // 598 //
833 G4FieldTrack RevisedEndPointFT = Cur << 599 G4FieldTrack newEndPointFT=
834 recalculatedB = << 600 ReEstimateEndpoint( CurrentA_PointVelocity,
835 CheckAndReEstimateEndpoint( Curre << 601 CurrentB_PointVelocity,
836 Curre << 602 linDistSq, // to avoid recalculation
837 Revis << 603 curveDist );
838 error << 604 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
839 if( recalculatedB ) << 605 CurrentB_PointVelocity = newEndPointFT;
840 { << 606 if (depth==1)
841 CurrentB_PointVelocity = RevisedEn << 607 {
842 << 608 recalculatedEndPoint = true;
843 if ( depth == 1 ) << 609 IntersectedOrRecalculatedFT = newEndPointFT;
844 { << 610 // 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 } 611 }
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 } 612 }
880 << 613
881 G4ThreeVector Point_A = CurrentA_Point << 614 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
882 G4ThreeVector Point_B = CurrentB_Point << 615 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
883 G4ThreeVector PointGi; <<
884 GetNavigatorFor()->LocateGlobalPointWi 616 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
885 G4bool Intersects_AB = IntersectChord( << 617 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
886 << 618 fPreviousSafety,
887 << 619 fPreviousSftOrigin,stepLengthAB,
888 << 620 PointGe);
889 if( Intersects_AB ) 621 if( Intersects_AB )
890 { 622 {
891 last_AF_intersection = Intersects_AB 623 last_AF_intersection = Intersects_AB;
892 CurrentE_Point = PointGi; << 624 CurrentE_Point = PointGe;
893 << 625 }
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--; 626 depth--;
906 fin_section_depth[depth] = true; << 627 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 } 628 }
927 #endif <<
928 } // if(!found_aproximate_intersection) 629 } // if(!found_aproximate_intersection)
929 630
930 assert( validIntersectP || there_is_no_int <<
931 || found_approxima <<
932 <<
933 } while ( ( ! found_approximate_intersection 631 } while ( ( ! found_approximate_intersection )
934 && ( ! there_is_no_intersection ) 632 && ( ! there_is_no_intersection )
935 && ( substep_no <= fMaxSteps) ); / << 633 && ( substep_no <= max_substeps) ); // UNTIL found or failed
936 634
937 if( substep_no > max_no_seen ) 635 if( substep_no > max_no_seen )
938 { 636 {
939 max_no_seen = substep_no; 637 max_no_seen = substep_no;
940 #ifdef G4DEBUG_FIELD << 638 if( max_no_seen > warn_substeps )
941 if( max_no_seen > fWarnSteps ) <<
942 { 639 {
943 trigger_substepno_print = max_no_seen-20 640 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
944 } 641 }
945 #endif <<
946 } 642 }
947 643
948 if( !there_is_no_intersection && !found_appr << 644 if( ( substep_no >= max_substeps)
>> 645 && !there_is_no_intersection
>> 646 && !found_approximate_intersection )
949 { 647 {
950 if( substep_no >= fMaxSteps) << 648 G4cerr << "WARNING - G4MultiLevelLocator::EstimateIntersectionPoint()"
951 { << 649 << G4endl
952 // Since we cannot go further (yet), w << 650 << " Convergence is requiring too many substeps: "
953 << 651 << substep_no << G4endl;
954 recalculatedEndPoint = true; << 652 G4cerr << " Abandoning effort to intersect. " << G4endl;
955 IntersectedOrRecalculatedFT = CurrentA << 653 G4cerr << " Information on start & current step follows in cout."
956 found_approximate_intersection = false << 654 << G4endl;
957 << 655 G4cout << "WARNING - G4MultiLevelLocator::EstimateIntersectionPoint()"
958 std::ostringstream message; << 656 << G4endl
959 message << G4endl; << 657 << " Convergence is requiring too many substeps: "
960 message << "Convergence is requiring t << 658 << substep_no << G4endl;
961 << substep_no << " ( limit = << 659 G4cout << " Found intersection = "
962 << G4endl << 660 << found_approximate_intersection << G4endl
963 << " Abandoning effort to int << 661 << " Intersection exists = "
964 message << " Number of calls to MLL << 662 << !there_is_no_intersection << G4endl;
965 message << " iteration = " << substep << 663 G4cout << " Start and Endpoint of Requested Step:" << G4endl;
966 << 664 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
967 message.precision( 12 ); << 665 -1.0, NewSafety, 0);
968 G4double done_len = CurrentA_PointVelo << 666 G4cout << G4endl;
969 G4double full_len = CurveEndPointVeloc << 667 G4cout << " 'Bracketing' starting and endpoint of current Sub-Step"
970 message << " Undertaken only le << 668 << G4endl;
971 << " out of " << full_len << " << 669 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
972 << " Remaining length = << 670 -1.0, NewSafety, substep_no-1);
973 << 671 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
974 message << " Start and end-point o << 672 -1.0, NewSafety, substep_no);
975 printStatus( CurveStartPointVelocity, << 673 G4cout << G4endl;
976 -1.0, NewSafety, 0, m << 674
977 message << " Start and end-poin << 675 #ifdef FUTURE_CORRECTION
978 printStatus( CurrentA_PointVelocity, C << 676 // 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 677
1014 void G4MultiLevelLocator::ReportFieldValue( c << 678 if ( ! found_approximate_intersection )
1015 c << 679 {
1016 c << 680 recalculatedEndPoint = true;
1017 { << 681 // Return the further valid intersection point -- potentially A ??
1018 enum { maxNumFieldComp = 24 }; << 682 // JA/19 Jan 2006
>> 683 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
>> 684
>> 685 G4cout << "WARNING - G4MultiLevelLocator::EstimateIntersectionPoint()"
>> 686 << G4endl
>> 687 << " Did not convergence after " << substep_no
>> 688 << " substeps." << G4endl;
>> 689 G4cout << " The endpoint was adjused to pointA resulting"
>> 690 << G4endl
>> 691 << " from the last substep: " << CurrentA_PointVelocity
>> 692 << G4endl;
>> 693 }
>> 694 #endif
1019 695
1020 G4ThreeVector position = locationPV.GetPos << 696 oldprc = G4cout.precision( 10 );
1021 G4double startPoint[4] = { position.x(), p << 697 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
1022 locationPV.GetL << 698 G4double full_len = CurveEndPointVelocity.GetCurveLength();
1023 G4double FieldVec[maxNumFieldComp]; // 24 << 699 G4cout << "ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
1024 for (double & i : FieldVec) << 700 << G4endl
1025 { << 701 << " Undertaken only length: " << done_len
1026 i = 0.0; << 702 << " out of " << full_len << " required." << G4endl;
1027 } << 703 G4cout << " Remaining length = " << full_len - done_len << G4endl;
1028 equation->GetFieldValue( startPoint, Field << 704 G4cout.precision( oldprc );
1029 G4cout << " B-field value (" << nameLoc < << 705
1030 << FieldVec[0] << " " << FieldVec[1 << 706 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
1031 G4double Emag2= G4ThreeVector( FieldVec[3] << 707 "UnableToLocateIntersection", FatalException,
1032 FieldVec[4] << 708 "Too many substeps while trying to locate intersection.");
1033 FieldVec[5] << 709 }
1034 if( Emag2 > 0.0 ) << 710 else if( substep_no >= warn_substeps )
1035 { << 711 {
1036 G4cout << " Electric = " << FieldVec[3] << 712 oldprc = G4cout.precision( 10 );
1037 << FieldVec[4] << 713 G4cout << "WARNING - G4MultiLevelLocator::EstimateIntersectionPoint()"
1038 << FieldVec[5] << 714 << G4endl
1039 } << 715 << " Undertaken length: "
1040 return; << 716 << CurrentB_PointVelocity.GetCurveLength();
>> 717 G4cout << " - Needed: " << substep_no << " substeps." << G4endl
>> 718 << " Warning level = " << warn_substeps
>> 719 << " and maximum substeps = " << max_substeps << G4endl;
>> 720 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
>> 721 "DifficultyToLocateIntersection", JustWarning,
>> 722 "Many substeps while trying to locate intersection.");
>> 723 G4cout.precision( oldprc );
>> 724 }
>> 725 return !there_is_no_intersection; // Success or failure
1041 } 726 }
1042 727