| 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: G4VIntersectionLocator.cc 102290 2017-01-20 11:19:44Z gcosmo $
>> 27 //
26 // Class G4VIntersectionLocator implementation 28 // Class G4VIntersectionLocator implementation
27 // 29 //
28 // 27.10.08 - John Apostolakis, Tatiana Nikiti 30 // 27.10.08 - John Apostolakis, Tatiana Nikitina.
29 // ------------------------------------------- 31 // ---------------------------------------------------------------------------
30 32
31 #include <iomanip> 33 #include <iomanip>
32 #include <sstream> 34 #include <sstream>
33 35
34 #include "globals.hh" 36 #include "globals.hh"
35 #include "G4ios.hh" 37 #include "G4ios.hh"
36 #include "G4AutoDelete.hh" 38 #include "G4AutoDelete.hh"
37 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh"
38 #include "G4VIntersectionLocator.hh" 40 #include "G4VIntersectionLocator.hh"
39 #include "G4TouchableHandle.hh" <<
40 #include "G4GeometryTolerance.hh" 41 #include "G4GeometryTolerance.hh"
41 42
42 ////////////////////////////////////////////// 43 ///////////////////////////////////////////////////////////////////////////
43 // 44 //
44 // Constructor 45 // Constructor
45 // 46 //
46 G4VIntersectionLocator::G4VIntersectionLocator << 47 G4VIntersectionLocator:: G4VIntersectionLocator(G4Navigator *theNavigator)
47 : fiNavigator(theNavigator) << 48 : fVerboseLevel(0),
>> 49 fUseNormalCorrection(false),
>> 50 fCheckMode(false),
>> 51 fiNavigator(theNavigator),
>> 52 fiChordFinder(0), // Not set - overridden at each step
>> 53 fiEpsilonStep(-1.0), // Out of range - overridden at each step
>> 54 fiDeltaIntersection(-1.0), // Out of range - overridden at each step
>> 55 fiUseSafety(false), // Default - overridden at each step
>> 56 fpTouchable(0)
48 { 57 {
49 kCarTolerance = G4GeometryTolerance::GetInst 58 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
50 << 59 fHelpingNavigator = new G4Navigator();
51 if( fiNavigator->GetExternalNavigation() == <<
52 { <<
53 fHelpingNavigator = new G4Navigator(); <<
54 } <<
55 else // Must clone the navigator, together w <<
56 { <<
57 fHelpingNavigator = fiNavigator->Clone(); <<
58 } <<
59 } 60 }
60 61
61 ////////////////////////////////////////////// 62 ///////////////////////////////////////////////////////////////////////////
62 // 63 //
63 // Destructor. 64 // Destructor.
64 // 65 //
65 G4VIntersectionLocator::~G4VIntersectionLocato 66 G4VIntersectionLocator::~G4VIntersectionLocator()
66 { 67 {
67 delete fHelpingNavigator; 68 delete fHelpingNavigator;
68 delete fpTouchable; 69 delete fpTouchable;
69 } 70 }
70 71
71 ////////////////////////////////////////////// 72 ///////////////////////////////////////////////////////////////////////////
72 // 73 //
73 // Dump status of propagator to cout (old meth 74 // Dump status of propagator to cout (old method)
74 // 75 //
75 void 76 void
76 G4VIntersectionLocator::printStatus( const G4F << 77 G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
77 const G4F << 78 const G4FieldTrack& CurrentFT,
78 G4d << 79 G4double requestStep,
79 G4d << 80 G4double safety,
80 G4i << 81 G4int stepNo)
81 { 82 {
82 std::ostringstream os; 83 std::ostringstream os;
83 printStatus( StartFT,CurrentFT,requestStep,s 84 printStatus( StartFT,CurrentFT,requestStep,safety,stepNo,os,fVerboseLevel);
84 G4cout << os.str(); 85 G4cout << os.str();
85 } 86 }
86 87
87 ////////////////////////////////////////////// 88 ///////////////////////////////////////////////////////////////////////////
88 // 89 //
89 // Dumps status of propagator. 90 // Dumps status of propagator.
90 // 91 //
91 void 92 void
92 G4VIntersectionLocator::printStatus( const G4F << 93 G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
93 const G4F << 94 const G4FieldTrack& CurrentFT,
94 G4d << 95 G4double requestStep,
95 G4d << 96 G4double safety,
96 G4i << 97 G4int stepNo,
97 std << 98 std::ostream& os,
98 G4i << 99 int verboseLevel)
99 { 100 {
100 // const G4int verboseLevel= fVerboseLevel; 101 // const G4int verboseLevel= fVerboseLevel;
101 const G4ThreeVector StartPosition = St 102 const G4ThreeVector StartPosition = StartFT.GetPosition();
102 const G4ThreeVector StartUnitVelocity = St 103 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
103 const G4ThreeVector CurrentPosition = Cu 104 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
104 const G4ThreeVector CurrentUnitVelocity = Cu 105 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
105 106
106 G4double step_len = CurrentFT.GetCurveLength 107 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
107 G4long oldprc; // cout/cerr precision setti << 108 G4int oldprc; // cout/cerr precision settings
108 109
109 if( ((stepNo == 0) && (verboseLevel <3)) || 110 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
110 { 111 {
111 oldprc = os.precision(4); 112 oldprc = os.precision(4);
112 os << std::setw( 6) << " " 113 os << std::setw( 6) << " "
113 << std::setw( 25) << " Current Posi 114 << std::setw( 25) << " Current Position and Direction" << " "
114 << G4endl; 115 << G4endl;
115 os << std::setw( 5) << "Step#" 116 os << std::setw( 5) << "Step#"
116 << std::setw(10) << " s " << " " 117 << std::setw(10) << " s " << " "
117 << std::setw(10) << "X(mm)" << " " 118 << std::setw(10) << "X(mm)" << " "
118 << std::setw(10) << "Y(mm)" << " " 119 << std::setw(10) << "Y(mm)" << " "
119 << std::setw(10) << "Z(mm)" << " " 120 << std::setw(10) << "Z(mm)" << " "
120 << std::setw( 7) << " N_x " << " " 121 << std::setw( 7) << " N_x " << " "
121 << std::setw( 7) << " N_y " << " " 122 << std::setw( 7) << " N_y " << " "
122 << std::setw( 7) << " N_z " << " " 123 << std::setw( 7) << " N_z " << " " ;
123 os << std::setw( 7) << " Delta|N|" << " " 124 os << std::setw( 7) << " Delta|N|" << " "
124 << std::setw( 9) << "StepLen" << " 125 << std::setw( 9) << "StepLen" << " "
125 << std::setw(12) << "StartSafety" < 126 << std::setw(12) << "StartSafety" << " "
126 << std::setw( 9) << "PhsStep" << " 127 << std::setw( 9) << "PhsStep" << " ";
127 os << G4endl; 128 os << G4endl;
128 os.precision(oldprc); 129 os.precision(oldprc);
129 } 130 }
130 if((stepNo == 0) && (verboseLevel <=3)) 131 if((stepNo == 0) && (verboseLevel <=3))
131 { 132 {
132 // Recurse to print the start values 133 // Recurse to print the start values
133 // 134 //
134 printStatus( StartFT, StartFT, -1.0, safet 135 printStatus( StartFT, StartFT, -1.0, safety, -1, os, verboseLevel);
135 } 136 }
136 if( verboseLevel <= 3 ) 137 if( verboseLevel <= 3 )
137 { 138 {
138 if( stepNo >= 0) 139 if( stepNo >= 0)
139 { 140 {
140 os << std::setw( 4) << stepNo << " "; 141 os << std::setw( 4) << stepNo << " ";
141 } 142 }
142 else 143 else
143 { 144 {
144 os << std::setw( 5) << "Start" ; 145 os << std::setw( 5) << "Start" ;
145 } 146 }
146 oldprc = os.precision(8); 147 oldprc = os.precision(8);
147 os << std::setw(10) << CurrentFT.GetCurveL 148 os << std::setw(10) << CurrentFT.GetCurveLength() << " ";
148 os << std::setw(10) << CurrentPosition.x() 149 os << std::setw(10) << CurrentPosition.x() << " "
149 << std::setw(10) << CurrentPosition 150 << std::setw(10) << CurrentPosition.y() << " "
150 << std::setw(10) << CurrentPosition 151 << std::setw(10) << CurrentPosition.z() << " ";
151 os.precision(4); 152 os.precision(4);
152 os << std::setw( 7) << CurrentUnitVelocity 153 os << std::setw( 7) << CurrentUnitVelocity.x() << " "
153 << std::setw( 7) << CurrentUnitVelo 154 << std::setw( 7) << CurrentUnitVelocity.y() << " "
154 << std::setw( 7) << CurrentUnitVelo 155 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
155 os.precision(3); 156 os.precision(3);
156 os << std::setw( 7) 157 os << std::setw( 7)
157 << CurrentFT.GetMomentum().mag()- S 158 << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag()
158 << " "; 159 << " ";
159 os << std::setw( 9) << step_len << " "; 160 os << std::setw( 9) << step_len << " ";
160 os << std::setw(12) << safety << " "; 161 os << std::setw(12) << safety << " ";
161 if( requestStep != -1.0 ) 162 if( requestStep != -1.0 )
162 { 163 {
163 os << std::setw( 9) << requestStep << " 164 os << std::setw( 9) << requestStep << " ";
164 } 165 }
165 else 166 else
166 { 167 {
167 os << std::setw( 9) << "Init/NotKnown" < 168 os << std::setw( 9) << "Init/NotKnown" << " ";
168 } 169 }
169 os << G4endl; 170 os << G4endl;
170 os.precision(oldprc); 171 os.precision(oldprc);
171 } 172 }
172 else // if( verboseLevel > 3 ) 173 else // if( verboseLevel > 3 )
173 { 174 {
174 // Multi-line output 175 // Multi-line output
175 176
176 os << "Step taken was " << step_len 177 os << "Step taken was " << step_len
177 << " out of PhysicalStep= " << req 178 << " out of PhysicalStep= " << requestStep << G4endl;
178 os << "Final safety is: " << safety << G4e 179 os << "Final safety is: " << safety << G4endl;
179 os << "Chord length = " << (CurrentPositio 180 os << "Chord length = " << (CurrentPosition-StartPosition).mag()
180 << G4endl; 181 << G4endl;
181 os << G4endl; 182 os << G4endl;
182 } 183 }
183 } 184 }
184 185
185 ////////////////////////////////////////////// 186 ///////////////////////////////////////////////////////////////////////////
186 // 187 //
187 // ReEstimateEndPoint. 188 // ReEstimateEndPoint.
188 // 189 //
189 G4FieldTrack G4VIntersectionLocator:: 190 G4FieldTrack G4VIntersectionLocator::
190 ReEstimateEndpoint( const G4FieldTrack& Curren 191 ReEstimateEndpoint( const G4FieldTrack& CurrentStateA,
191 const G4FieldTrack& Estima 192 const G4FieldTrack& EstimatedEndStateB,
>> 193 // bool & recalculated)
192 G4double , // l 194 G4double , // linearDistSq, // NOT used
193 G4double << 195 G4double ) // curveDist ) // NOT used
194 #ifdef G4DEBUG_FIELD <<
195 curveDist <<
196 #endif <<
197 ) <<
198 { 196 {
199 G4FieldTrack newEndPoint( CurrentStateA ); 197 G4FieldTrack newEndPoint( CurrentStateA );
200 auto integrDriver = GetChordFinderFor()->Get << 198 G4MagInt_Driver* integrDriver = GetChordFinderFor()->GetIntegrationDriver();
201 199
202 G4FieldTrack retEndPoint( CurrentStateA ); 200 G4FieldTrack retEndPoint( CurrentStateA );
203 G4bool goodAdvance; 201 G4bool goodAdvance;
204 G4int itrial = 0; << 202 G4int itrial=0;
205 const G4int no_trials = 20; << 203 const G4int no_trials=20;
206 204
207 205
208 G4double endCurveLen= EstimatedEndStateB.Get 206 G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
209 207
210 do // Loop checking, 07.10.2016, JA << 208 do // Loop checking, 07.10.2016, J.Apostolakis
211 { 209 {
212 G4double currentCurveLen = newEndPoint.Get << 210 G4double currentCurveLen= newEndPoint.GetCurveLength();
213 G4double advanceLength = endCurveLen - cur << 211 G4double advanceLength= endCurveLen - currentCurveLen ;
214 if (std::abs(advanceLength)<kCarTolerance) 212 if (std::abs(advanceLength)<kCarTolerance)
215 { 213 {
216 goodAdvance=true; 214 goodAdvance=true;
217 } 215 }
218 else 216 else
219 { 217 {
220 goodAdvance = integrDriver->AccurateAdva << 218 goodAdvance= integrDriver->AccurateAdvance(newEndPoint, advanceLength,
221 << 219 GetEpsilonStepFor());
222 } << 220 }
223 } 221 }
224 while( !goodAdvance && (++itrial < no_trials 222 while( !goodAdvance && (++itrial < no_trials) );
225 223
226 if( goodAdvance ) 224 if( goodAdvance )
227 { 225 {
228 retEndPoint = newEndPoint; 226 retEndPoint = newEndPoint;
229 } 227 }
230 else 228 else
231 { 229 {
232 retEndPoint = EstimatedEndStateB; // Could 230 retEndPoint = EstimatedEndStateB; // Could not improve without major work !!
233 } 231 }
234 232
235 // All the work is done 233 // All the work is done
236 // below are some diagnostics only -- befor 234 // below are some diagnostics only -- before the return!
237 // 235 //
238 const G4String MethodName("G4VIntersectionLo 236 const G4String MethodName("G4VIntersectionLocator::ReEstimateEndpoint()");
239 237
240 #ifdef G4VERBOSE 238 #ifdef G4VERBOSE
241 G4int latest_good_trials = 0; 239 G4int latest_good_trials = 0;
242 if( itrial > 1) 240 if( itrial > 1)
243 { 241 {
244 if( fVerboseLevel > 0 ) 242 if( fVerboseLevel > 0 )
245 { 243 {
246 G4cout << MethodName << " called - goodA 244 G4cout << MethodName << " called - goodAdv= " << goodAdvance
247 << " trials = " << itrial 245 << " trials = " << itrial
248 << " previous good= " << latest_g 246 << " previous good= " << latest_good_trials
249 << G4endl; 247 << G4endl;
250 } 248 }
251 latest_good_trials = 0; 249 latest_good_trials = 0;
252 } 250 }
253 else 251 else
254 { 252 {
255 ++latest_good_trials; << 253 latest_good_trials++;
256 } 254 }
257 #endif 255 #endif
258 256
259 #ifdef G4DEBUG_FIELD 257 #ifdef G4DEBUG_FIELD
260 G4double lengthDone = newEndPoint.GetCurveLe 258 G4double lengthDone = newEndPoint.GetCurveLength()
261 - CurrentStateA.GetCurve 259 - CurrentStateA.GetCurveLength();
262 if( !goodAdvance ) 260 if( !goodAdvance )
263 { 261 {
264 if( fVerboseLevel >= 3 ) 262 if( fVerboseLevel >= 3 )
265 { 263 {
266 G4cout << MethodName << "> AccurateAdvan 264 G4cout << MethodName << "> AccurateAdvance failed " ;
267 G4cout << " in " << itrial << " integrat 265 G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
268 G4cout << " It went only " << lengthDone 266 G4cout << " It went only " << lengthDone << " instead of " << curveDist
269 << " -- a difference of " << curv 267 << " -- a difference of " << curveDist - lengthDone << G4endl;
270 G4cout << " ReEstimateEndpoint> Reset en 268 G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
271 << G4endl; 269 << G4endl;
272 } 270 }
273 } 271 }
274 G4double linearDist = ( EstimatedEndStateB.G 272 G4double linearDist = ( EstimatedEndStateB.GetPosition()
275 - CurrentStateA.GetPosit 273 - CurrentStateA.GetPosition() ).mag();
276 static G4int noInaccuracyWarnings = 0; 274 static G4int noInaccuracyWarnings = 0;
277 G4int maxNoWarnings = 10; 275 G4int maxNoWarnings = 10;
278 if ( (noInaccuracyWarnings < maxNoWarnings 276 if ( (noInaccuracyWarnings < maxNoWarnings )
279 || (fVerboseLevel > 1) ) 277 || (fVerboseLevel > 1) )
280 { 278 {
281 G4ThreeVector move = newEndPoint.GetPosi 279 G4ThreeVector move = newEndPoint.GetPosition()
282 - EstimatedEndStateB. 280 - EstimatedEndStateB.GetPosition();
283 std::ostringstream message; 281 std::ostringstream message;
284 message.precision(12); 282 message.precision(12);
285 message << " Integration inaccuracy requ 283 message << " Integration inaccuracy requires"
286 << " an adjustment in the step's 284 << " an adjustment in the step's endpoint." << G4endl
287 << " Two mid-points are furthe 285 << " Two mid-points are further apart than their"
288 << " curve length difference" 286 << " curve length difference" << G4endl
289 << " Dist = " << linearD 287 << " Dist = " << linearDist
290 << " curve length = " << curveDi 288 << " curve length = " << curveDist << G4endl;
291 message << " Correction applied is " << 289 message << " Correction applied is " << move.mag() << G4endl
292 << " Old Estimated B position= 290 << " Old Estimated B position= "
293 << EstimatedEndStateB.GetPositio 291 << EstimatedEndStateB.GetPosition() << G4endl
294 << " Recalculated Position= 292 << " Recalculated Position= "
295 << newEndPoint.GetPosition() << 293 << newEndPoint.GetPosition() << G4endl
296 << " Change ( new - old ) = 294 << " Change ( new - old ) = " << move;
297 G4Exception("G4VIntersectionLocator::ReE 295 G4Exception("G4VIntersectionLocator::ReEstimateEndpoint()",
298 "GeomNav1002", JustWarning, 296 "GeomNav1002", JustWarning, message);
299 } 297 }
300 /* <<
301 #else 298 #else
302 // Statistics on the RMS value of the correc 299 // Statistics on the RMS value of the corrections
303 300
304 static G4ThreadLocal G4int noCorrections = 0 << 301 static G4ThreadLocal G4int noCorrections=0;
305 ++noCorrections; << 302 static G4ThreadLocal G4double sumCorrectionsSq = 0;
>> 303 noCorrections++;
306 if( goodAdvance ) 304 if( goodAdvance )
307 { 305 {
308 static G4ThreadLocal G4double sumCorrectio <<
309 sumCorrectionsSq += (EstimatedEndStateB.Ge 306 sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
310 newEndPoint.GetPositi 307 newEndPoint.GetPosition()).mag2();
311 } 308 }
312 */ <<
313 #endif 309 #endif
314 310
315 return retEndPoint; 311 return retEndPoint;
316 } 312 }
317 313
318 314
319 ////////////////////////////////////////////// 315 ///////////////////////////////////////////////////////////////////////////
320 // 316 //
321 // ReEstimateEndPoint. 317 // ReEstimateEndPoint.
322 // 318 //
323 // The following values are returned in curv 319 // The following values are returned in curveError
324 // 0: Normal - no problem 320 // 0: Normal - no problem
325 // 1: Unexpected co-incidence - milder m 321 // 1: Unexpected co-incidence - milder mixup
326 // 2: Real mixup - EndB is NOT past Star 322 // 2: Real mixup - EndB is NOT past StartA
327 // ( ie. StartA's curve-lengh is bi 323 // ( ie. StartA's curve-lengh is bigger than EndB's)
328 324
329 325
330 G4bool G4VIntersectionLocator:: 326 G4bool G4VIntersectionLocator::
331 CheckAndReEstimateEndpoint( const G4FieldTrack 327 CheckAndReEstimateEndpoint( const G4FieldTrack& CurrentStartA,
332 const G4FieldTrack 328 const G4FieldTrack& EstimatedEndB,
333 G4FieldTrack 329 G4FieldTrack& RevisedEndPoint,
334 G4int& << 330 G4int & curveError)
335 { 331 {
336 G4double linDistSq, curveDist; 332 G4double linDistSq, curveDist;
337 333
338 G4bool recalculated = false; << 334 G4bool recalculated= false;
339 curveError= 0; 335 curveError= 0;
340 336
341 linDistSq = ( EstimatedEndB.GetPosition() 337 linDistSq = ( EstimatedEndB.GetPosition()
342 - CurrentStartA.GetPosition() ). << 338 - CurrentStartA.GetPosition() ).mag2();
343 curveDist = EstimatedEndB.GetCurveLength() 339 curveDist = EstimatedEndB.GetCurveLength()
344 - CurrentStartA.GetCurveLength(); << 340 - CurrentStartA.GetCurveLength();
345 if( (curveDist>=0.0) 341 if( (curveDist>=0.0)
346 && (curveDist*curveDist *(1.0+2.0*fiEpsil 342 && (curveDist*curveDist *(1.0+2.0*fiEpsilonStep ) < linDistSq ) )
347 { 343 {
348 G4FieldTrack newEndPointFT = EstimatedEndB << 344 // G4FieldTrack oldPointVelB = EstimatedEndB;
>> 345 G4FieldTrack newEndPointFT= EstimatedEndB; // Unused
349 346
350 if (curveDist>0.0) 347 if (curveDist>0.0)
351 { 348 {
352 // Re-integrate to obtain a new B 349 // Re-integrate to obtain a new B
353 RevisedEndPoint = ReEstimateEndpoint( C << 350 RevisedEndPoint= ReEstimateEndpoint( CurrentStartA,
354 E << 351 EstimatedEndB,
355 l << 352 linDistSq,
356 c << 353 curveDist );
357 recalculated = true; 354 recalculated = true;
358 } 355 }
359 else 356 else
360 { 357 {
361 // Zero length -> no advance! 358 // Zero length -> no advance!
362 newEndPointFT = CurrentStartA; << 359 newEndPointFT= CurrentStartA;
363 recalculated = true; 360 recalculated = true;
364 curveError = 1; // Unexpected co-incid 361 curveError = 1; // Unexpected co-incidence - milder mixup
365 362
366 G4Exception("G4MultiLevelLocator::Estim 363 G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
367 "GeomNav1002", JustWarning, 364 "GeomNav1002", JustWarning,
368 "A & B are at equal distance in 2nd 365 "A & B are at equal distance in 2nd half. A & B will coincide." );
369 } 366 }
370 } 367 }
371 368
372 // Sanity check 369 // Sanity check
373 // 370 //
374 if( curveDist < 0.0 ) 371 if( curveDist < 0.0 )
375 { 372 {
>> 373 // clean = false;
376 curveError = 2; // Real mixup 374 curveError = 2; // Real mixup
377 } 375 }
378 return recalculated; 376 return recalculated;
379 } 377 }
380 378
381 ////////////////////////////////////////////// 379 ///////////////////////////////////////////////////////////////////////////
382 // 380 //
383 // Method for finding SurfaceNormal of Interse 381 // Method for finding SurfaceNormal of Intersecting Solid
384 // 382 //
385 G4ThreeVector G4VIntersectionLocator:: 383 G4ThreeVector G4VIntersectionLocator::
386 GetLocalSurfaceNormal(const G4ThreeVector& Cur 384 GetLocalSurfaceNormal(const G4ThreeVector& CurrentE_Point, G4bool& validNormal)
387 { 385 {
388 G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0 386 G4ThreeVector Normal(G4ThreeVector(0.0,0.0,0.0));
389 G4VPhysicalVolume* located; 387 G4VPhysicalVolume* located;
390 388
391 validNormal = false; 389 validNormal = false;
392 fHelpingNavigator->SetWorldVolume(GetNavigat 390 fHelpingNavigator->SetWorldVolume(GetNavigatorFor()->GetWorldVolume());
393 located = fHelpingNavigator->LocateGlobalPoi 391 located = fHelpingNavigator->LocateGlobalPointAndSetup( CurrentE_Point );
394 392
395 delete fpTouchable; 393 delete fpTouchable;
396 fpTouchable = fHelpingNavigator->CreateTouch 394 fpTouchable = fHelpingNavigator->CreateTouchableHistory();
397 395
398 // To check if we can use GetGlobalExitNorma 396 // To check if we can use GetGlobalExitNormal()
399 // 397 //
400 G4ThreeVector localPosition = fpTouchable->G 398 G4ThreeVector localPosition = fpTouchable->GetHistory()
401 ->GetTopTransform().TransformP 399 ->GetTopTransform().TransformPoint(CurrentE_Point);
402 400
403 // Issue: in the case of coincident surfaces 401 // Issue: in the case of coincident surfaces, this version does not recognise
404 // which side you are located onto (c 402 // which side you are located onto (can return vector with wrong sign.)
405 // TO-DO: use direction (of chord) to identi 403 // TO-DO: use direction (of chord) to identify volume we will be "entering"
406 404
407 if( located != nullptr) << 405 if( located != 0)
408 { 406 {
409 G4LogicalVolume* pLogical= located->GetLog 407 G4LogicalVolume* pLogical= located->GetLogicalVolume();
410 G4VSolid* pSolid; 408 G4VSolid* pSolid;
411 409
412 if( (pLogical != nullptr) && ( (pSolid=pLo << 410 if( (pLogical != 0) && ( (pSolid=pLogical->GetSolid()) !=0 ) )
413 { 411 {
>> 412 // G4bool goodPoint, nearbyPoint;
>> 413 // G4int numGoodPoints, numNearbyPoints; // --> use for stats
414 if ( ( pSolid->Inside(localPosition)==kS 414 if ( ( pSolid->Inside(localPosition)==kSurface )
415 || ( pSolid->DistanceToOut(localPositi << 415 || ( pSolid->DistanceToOut(localPosition) < 1000.0 * kCarTolerance )
>> 416 )
416 { 417 {
417 Normal = pSolid->SurfaceNormal(localPo 418 Normal = pSolid->SurfaceNormal(localPosition);
418 validNormal = true; 419 validNormal = true;
419 420
420 #ifdef G4DEBUG_FIELD 421 #ifdef G4DEBUG_FIELD
421 if( std::fabs(Normal.mag2() - 1.0 ) > 422 if( std::fabs(Normal.mag2() - 1.0 ) > CLHEP::perThousand)
422 { 423 {
423 G4cerr << "PROBLEM in G4VIntersectio 424 G4cerr << "PROBLEM in G4VIntersectionLocator::GetLocalSurfaceNormal."
424 << G4endl; 425 << G4endl;
425 G4cerr << " Normal is not unit - ma 426 G4cerr << " Normal is not unit - mag=" << Normal.mag() << G4endl;
426 G4cerr << " at trial local point " 427 G4cerr << " at trial local point " << CurrentE_Point << G4endl;
427 G4cerr << " Solid is " << *pSolid 428 G4cerr << " Solid is " << *pSolid << G4endl;
428 } 429 }
429 #endif 430 #endif
430 } 431 }
431 } 432 }
432 } 433 }
433 return Normal; 434 return Normal;
434 } 435 }
435 436
436 ////////////////////////////////////////////// 437 ///////////////////////////////////////////////////////////////////////////
437 // 438 //
438 // Adjustment of Found Intersection 439 // Adjustment of Found Intersection
439 // 440 //
440 G4bool G4VIntersectionLocator:: 441 G4bool G4VIntersectionLocator::
441 AdjustmentOfFoundIntersection( const G4ThreeVe 442 AdjustmentOfFoundIntersection( const G4ThreeVector& CurrentA_Point,
442 const G4ThreeVe 443 const G4ThreeVector& CurrentE_Point,
443 const G4ThreeVe 444 const G4ThreeVector& CurrentF_Point,
444 const G4ThreeVe 445 const G4ThreeVector& MomentumDir,
445 const G4bool 446 const G4bool IntersectAF,
446 G4ThreeVe 447 G4ThreeVector& IntersectionPoint, // I/O
447 G4double& 448 G4double& NewSafety, // I/O
448 G4double& 449 G4double& fPreviousSafety, // I/O
449 G4ThreeVe 450 G4ThreeVector& fPreviousSftOrigin )// I/O
450 { 451 {
451 G4double dist,lambda; 452 G4double dist,lambda;
452 G4ThreeVector Normal, NewPoint, Point_G; 453 G4ThreeVector Normal, NewPoint, Point_G;
453 G4bool goodAdjust = false, Intersects_FP = f << 454 G4bool goodAdjust=false, Intersects_FP=false, validNormal=false;
454 455
455 // Get SurfaceNormal of Intersecting Solid 456 // Get SurfaceNormal of Intersecting Solid
456 // 457 //
457 Normal = GetGlobalSurfaceNormal(CurrentE_Poi 458 Normal = GetGlobalSurfaceNormal(CurrentE_Point,validNormal);
458 if(!validNormal) { return false; } 459 if(!validNormal) { return false; }
459 460
460 // Intersection between Line and Plane 461 // Intersection between Line and Plane
461 // 462 //
462 G4double n_d_m = Normal.dot(MomentumDir); 463 G4double n_d_m = Normal.dot(MomentumDir);
463 if ( std::abs(n_d_m)>kCarTolerance ) 464 if ( std::abs(n_d_m)>kCarTolerance )
464 { 465 {
465 #ifdef G4VERBOSE 466 #ifdef G4VERBOSE
466 if ( fVerboseLevel>1 ) 467 if ( fVerboseLevel>1 )
467 { 468 {
468 G4Exception("G4VIntersectionLocator::Adj << 469 G4cerr << "WARNING - "
469 "GeomNav0003", JustWarning, << 470 << "G4VIntersectionLocator::AdjustementOfFoundIntersection()"
470 "No intersection. Parallels << 471 << G4endl
>> 472 << " No intersection. Parallels lines!" << G4endl;
471 } 473 }
472 #endif 474 #endif
473 lambda =- Normal.dot(CurrentF_Point-Curren 475 lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
474 476
475 // New candidate for Intersection 477 // New candidate for Intersection
476 // 478 //
477 NewPoint = CurrentF_Point+lambda*MomentumD 479 NewPoint = CurrentF_Point+lambda*MomentumDir;
478 480
479 // Distance from CurrentF to Calculated In 481 // Distance from CurrentF to Calculated Intersection
480 // 482 //
481 dist = std::abs(lambda); 483 dist = std::abs(lambda);
482 484
483 if ( dist<kCarTolerance*0.001 ) { return 485 if ( dist<kCarTolerance*0.001 ) { return false; }
484 486
485 // Calculation of new intersection point o 487 // Calculation of new intersection point on the path.
486 // 488 //
487 if ( IntersectAF ) // First part interse 489 if ( IntersectAF ) // First part intersects
488 { 490 {
489 G4double stepLengthFP; 491 G4double stepLengthFP;
490 G4ThreeVector Point_P = CurrentA_Point; 492 G4ThreeVector Point_P = CurrentA_Point;
491 GetNavigatorFor()->LocateGlobalPointWith 493 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_P);
492 Intersects_FP = IntersectChord( Point_P, 494 Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
493 fPreviou 495 fPreviousSafety, fPreviousSftOrigin,
494 stepLeng 496 stepLengthFP, Point_G );
495 497
496 } 498 }
497 else // Second part intersects 499 else // Second part intersects
498 { 500 {
499 G4double stepLengthFP; 501 G4double stepLengthFP;
500 GetNavigatorFor()->LocateGlobalPointWith 502 GetNavigatorFor()->LocateGlobalPointWithinVolume(CurrentF_Point );
501 Intersects_FP = IntersectChord( CurrentF 503 Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
502 fPreviou 504 fPreviousSafety, fPreviousSftOrigin,
503 stepLeng 505 stepLengthFP, Point_G );
504 } 506 }
505 if ( Intersects_FP ) 507 if ( Intersects_FP )
506 { 508 {
507 goodAdjust = true; 509 goodAdjust = true;
508 IntersectionPoint = Point_G; 510 IntersectionPoint = Point_G;
509 } 511 }
510 } 512 }
511 513
512 return goodAdjust; 514 return goodAdjust;
513 } 515 }
514 516
515 ////////////////////////////////////////////// 517 ///////////////////////////////////////////////////////////////////////////
516 // 518 //
517 // GetSurfaceNormal. 519 // GetSurfaceNormal.
518 // 520 //
519 G4ThreeVector G4VIntersectionLocator:: 521 G4ThreeVector G4VIntersectionLocator::
520 GetSurfaceNormal(const G4ThreeVector& CurrentI 522 GetSurfaceNormal(const G4ThreeVector& CurrentInt_Point,
521 G4bool& validNormal) << 523 G4bool& validNormal) // const
522 { 524 {
523 G4ThreeVector NormalAtEntry; // ( -10. , -10 525 G4ThreeVector NormalAtEntry; // ( -10. , -10., -10. );
524 526
525 G4ThreeVector NormalAtEntryLast, NormalAtEnt 527 G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
526 G4bool validNormalLast; 528 G4bool validNormalLast;
527 529
528 // Relies on a call to Navigator::ComputeSte 530 // Relies on a call to Navigator::ComputeStep in IntersectChord before
529 // this call 531 // this call
530 // 532 //
531 NormalAtEntryLast = GetLastSurfaceNormal( Cu 533 NormalAtEntryLast = GetLastSurfaceNormal( CurrentInt_Point, validNormalLast );
532 // May return valid=false in cases, includ 534 // May return valid=false in cases, including
533 // - if the candidate volume was not foun 535 // - if the candidate volume was not found (eg exiting world), or
534 // - a replica was involved -- determined 536 // - a replica was involved -- determined the step size.
535 // (This list is not complete.) 537 // (This list is not complete.)
536 538
537 #ifdef G4DEBUG_FIELD 539 #ifdef G4DEBUG_FIELD
538 if ( validNormalLast 540 if ( validNormalLast
539 && ( std::fabs(NormalAtEntryLast.mag2() - 1 541 && ( std::fabs(NormalAtEntryLast.mag2() - 1.0) > perThousand ) )
540 { 542 {
541 std::ostringstream message; 543 std::ostringstream message;
542 message << "PROBLEM: Normal is not unit - 544 message << "PROBLEM: Normal is not unit - magnitude = "
543 << NormalAtEntryLast.mag() << G4en 545 << NormalAtEntryLast.mag() << G4endl;
544 message << " at trial intersection point 546 message << " at trial intersection point " << CurrentInt_Point << G4endl;
545 message << " Obtained from Get *Last* Su 547 message << " Obtained from Get *Last* Surface Normal.";
546 G4Exception("G4VIntersectionLocator::GetSu 548 G4Exception("G4VIntersectionLocator::GetSurfaceNormal()",
547 "GeomNav1002", JustWarning, me 549 "GeomNav1002", JustWarning, message);
548 } 550 }
549 #endif 551 #endif
550 552
551 if( validNormalLast ) 553 if( validNormalLast )
552 { 554 {
553 NormalAtEntry = NormalAtEntryLast; << 555 NormalAtEntry=NormalAtEntryLast;
554 } 556 }
555 validNormal = validNormalLast; << 557 validNormal = validNormalLast;
556 558
557 return NormalAtEntry; 559 return NormalAtEntry;
558 } 560 }
559 561
560 ////////////////////////////////////////////// 562 ///////////////////////////////////////////////////////////////////////////
561 // 563 //
562 // GetGlobalSurfaceNormal. 564 // GetGlobalSurfaceNormal.
563 // 565 //
564 G4ThreeVector G4VIntersectionLocator:: 566 G4ThreeVector G4VIntersectionLocator::
565 GetGlobalSurfaceNormal(const G4ThreeVector& Cu 567 GetGlobalSurfaceNormal(const G4ThreeVector& CurrentE_Point,
566 G4bool& validNorm 568 G4bool& validNormal)
567 { 569 {
568 G4ThreeVector localNormal = GetLocalSurfaceN << 570 G4ThreeVector localNormal=
569 G4AffineTransform localToGlobal = / << 571 GetLocalSurfaceNormal( CurrentE_Point, validNormal );
>> 572 G4AffineTransform localToGlobal= // Must use the same Navigator !!
570 fHelpingNavigator->GetLocalToGlobalTrans 573 fHelpingNavigator->GetLocalToGlobalTransform();
571 G4ThreeVector globalNormal = localToGlobal.T << 574 G4ThreeVector globalNormal =
>> 575 localToGlobal.TransformAxis( localNormal );
572 576
573 #ifdef G4DEBUG_FIELD 577 #ifdef G4DEBUG_FIELD
574 if( validNormal && ( std::fabs(globalNormal. 578 if( validNormal && ( std::fabs(globalNormal.mag2() - 1.0) > perThousand ) )
575 { 579 {
576 std::ostringstream message; 580 std::ostringstream message;
577 message << "****************************** 581 message << "**************************************************************"
578 << G4endl; 582 << G4endl;
579 message << " Bad Normal in G4VIntersection 583 message << " Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
580 << G4endl; 584 << G4endl;
581 message << " * Constituents: " << G4endl; 585 message << " * Constituents: " << G4endl;
582 message << " Local Normal= " << localN 586 message << " Local Normal= " << localNormal << G4endl;
583 message << " Transform: " << G4endl 587 message << " Transform: " << G4endl
584 << " Net Translation= " << lo 588 << " Net Translation= " << localToGlobal.NetTranslation()
585 << G4endl 589 << G4endl
586 << " Net Rotation = " << lo 590 << " Net Rotation = " << localToGlobal.NetRotation()
587 << G4endl; 591 << G4endl;
588 message << " * Result: " << G4endl; 592 message << " * Result: " << G4endl;
589 message << " Global Normal= " << local 593 message << " Global Normal= " << localNormal << G4endl;
590 message << "****************************** << 594 message << "**************************************************************"
>> 595 << G4endl;
591 G4Exception("G4VIntersectionLocator::GetGl 596 G4Exception("G4VIntersectionLocator::GetGlobalSurfaceNormal()",
592 "GeomNav1002", JustWarning, me 597 "GeomNav1002", JustWarning, message);
593 } 598 }
594 #endif 599 #endif
595 600
596 return globalNormal; 601 return globalNormal;
597 } 602 }
598 603
599 ////////////////////////////////////////////// 604 ///////////////////////////////////////////////////////////////////////////
600 // 605 //
601 // GetLastSurfaceNormal. 606 // GetLastSurfaceNormal.
602 // 607 //
603 G4ThreeVector G4VIntersectionLocator:: 608 G4ThreeVector G4VIntersectionLocator::
604 GetLastSurfaceNormal( const G4ThreeVector& int 609 GetLastSurfaceNormal( const G4ThreeVector& intersectPoint,
605 G4bool& normalIsVa 610 G4bool& normalIsValid) const
606 { 611 {
607 G4ThreeVector normalVec; 612 G4ThreeVector normalVec;
608 G4bool validNorm; << 613 G4bool validNorm;
609 normalVec = fiNavigator->GetGlobalExitNormal 614 normalVec = fiNavigator->GetGlobalExitNormal( intersectPoint, &validNorm );
610 normalIsValid = validNorm; << 615 normalIsValid= validNorm;
611 616
612 return normalVec; 617 return normalVec;
613 } 618 }
614 619
615 ////////////////////////////////////////////// 620 ///////////////////////////////////////////////////////////////////////////
616 // 621 //
617 // ReportTrialStep. 622 // ReportTrialStep.
618 // 623 //
619 void G4VIntersectionLocator::ReportTrialStep( 624 void G4VIntersectionLocator::ReportTrialStep( G4int step_no,
620 const 625 const G4ThreeVector& ChordAB_v,
621 const 626 const G4ThreeVector& ChordEF_v,
622 const 627 const G4ThreeVector& NewMomentumDir,
623 const 628 const G4ThreeVector& NormalAtEntry,
624 629 G4bool validNormal )
625 { 630 {
626 G4double ABchord_length = ChordAB_v.mag(); << 631 G4double ABchord_length = ChordAB_v.mag();
627 G4double MomDir_dot_Norm = NewMomentumDir.do << 632 G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ) ;
628 G4double MomDir_dot_ABchord; << 633 G4double MomDir_dot_ABchord;
629 MomDir_dot_ABchord = (1.0 / ABchord_length) << 634 MomDir_dot_ABchord= (1.0 / ABchord_length) * NewMomentumDir.dot( ChordAB_v );
630 635
631 std::ostringstream outStream; 636 std::ostringstream outStream;
632 outStream << std::setw(6) << " Step# " << 637 outStream // G4cout
>> 638 << std::setw(6) << " Step# "
633 << std::setw(17) << " |ChordEF|(mag)" << " 639 << std::setw(17) << " |ChordEF|(mag)" << " "
634 << std::setw(18) << " uMomentum.Normal" << 640 << std::setw(18) << " uMomentum.Normal" << " "
635 << std::setw(18) << " uMomentum.ABdir " << 641 << std::setw(18) << " uMomentum.ABdir " << " "
636 << std::setw(16) << " AB-dist " << 642 << std::setw(16) << " AB-dist " << " "
637 << " Chord Vector (EF) " 643 << " Chord Vector (EF) "
638 << G4endl; 644 << G4endl;
639 outStream.precision(7); 645 outStream.precision(7);
640 outStream << " " << std::setw(5) << step_no << 646 outStream // G4cout
>> 647 << " " << std::setw(5) << step_no
641 << " " << std::setw(18) << ChordEF_v.mag() 648 << " " << std::setw(18) << ChordEF_v.mag()
642 << " " << std::setw(18) << MomDir_dot_Norm 649 << " " << std::setw(18) << MomDir_dot_Norm
643 << " " << std::setw(18) << MomDir_dot_ABch 650 << " " << std::setw(18) << MomDir_dot_ABchord
644 << " " << std::setw(12) << ABchord_length 651 << " " << std::setw(12) << ABchord_length
645 << " " << ChordEF_v 652 << " " << ChordEF_v
646 << G4endl; 653 << G4endl;
647 outStream << " MomentumDir= " << " " << NewM << 654 outStream // G4cout
>> 655 << " MomentumDir= " << " " << NewMomentumDir
648 << " Normal at Entry E= " << NormalAtEntry 656 << " Normal at Entry E= " << NormalAtEntry
649 << " AB chord = " << ChordAB_v 657 << " AB chord = " << ChordAB_v
650 << G4endl; 658 << G4endl;
651 G4cout << outStream.str(); << 659 G4cout << outStream.str(); // ostr_verbose;
652 660
653 if( ( std::fabs(NormalAtEntry.mag2() - 1.0) 661 if( ( std::fabs(NormalAtEntry.mag2() - 1.0) > perThousand ) )
654 { 662 {
655 std::ostringstream message; << 663 G4cerr << " PROBLEM in G4VIntersectionLocator::ReportTrialStep " << G4endl
656 message << "Normal is not unit - mag= " << << 664 << " Normal is not unit - mag=" << NormalAtEntry.mag()
657 << " ValidNormalAtE = " << << 665 << " ValidNormalAtE = " << validNormal
658 G4Exception("G4VIntersectionLocator::Repor << 666 << G4endl;
659 "GeomNav1002", JustWarning, me <<
660 } 667 }
661 return; 668 return;
662 } 669 }
663 670
664 ////////////////////////////////////////////// 671 ///////////////////////////////////////////////////////////////////////////
665 // 672 //
666 // LocateGlobalPointWithinVolumeAndCheck. 673 // LocateGlobalPointWithinVolumeAndCheck.
667 // 674 //
668 // Locate point using navigator: updates state 675 // Locate point using navigator: updates state of Navigator
669 // By default, it assumes that the point is in 676 // By default, it assumes that the point is inside the current volume,
670 // and returns true. 677 // and returns true.
671 // In check mode, checks whether the point is 678 // In check mode, checks whether the point is *inside* the volume.
672 // If it is inside, it returns true 679 // If it is inside, it returns true
673 // If not, issues a warning and returns fals 680 // If not, issues a warning and returns false.
674 // 681 //
675 G4bool G4VIntersectionLocator:: 682 G4bool G4VIntersectionLocator::
676 LocateGlobalPointWithinVolumeAndCheck( const G 683 LocateGlobalPointWithinVolumeAndCheck( const G4ThreeVector& position )
677 { 684 {
678 G4bool good = true; << 685 G4bool good= true;
679 G4Navigator* nav = GetNavigatorFor(); << 686 G4Navigator* nav= GetNavigatorFor();
680 const G4String 687 const G4String
681 MethodName("G4VIntersectionLocator::LocateGl 688 MethodName("G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()");
682 689
683 if( fCheckMode ) 690 if( fCheckMode )
684 { 691 {
685 G4bool navCheck= nav->IsCheckModeActive(); 692 G4bool navCheck= nav->IsCheckModeActive(); // Recover original value
686 nav->CheckMode(true); 693 nav->CheckMode(true);
687 694
688 // Identify the current volume 695 // Identify the current volume
689 696
690 G4TouchableHandle startTH= nav->CreateTouc << 697 G4TouchableHistoryHandle startTH= nav->CreateTouchableHistoryHandle();
691 G4VPhysicalVolume* motherPhys = startTH->G << 698 G4VPhysicalVolume* motherPhys= startTH->GetVolume();
692 G4VSolid* motherSolid = startTH-> << 699 G4VSolid* motherSolid= startTH->GetSolid();
693 G4AffineTransform transform = nav->GetGlob 700 G4AffineTransform transform = nav->GetGlobalToLocalTransform();
694 G4int motherCopyNo = motherPhys->GetCopyNo << 701 // GetLocalToGlobalTransform();
>> 702 G4int motherCopyNo= motherPhys->GetCopyNo();
695 703
696 // Let's check that the point is inside th 704 // Let's check that the point is inside the current solid
697 G4ThreeVector localPosition = transform.T 705 G4ThreeVector localPosition = transform.TransformPoint(position);
698 EInside inMother = motherSolid->Ins << 706 EInside inMother= motherSolid->Inside( localPosition );
699 if( inMother != kInside ) 707 if( inMother != kInside )
700 { 708 {
701 std::ostringstream message; << 709 G4cerr << " ERROR in " << MethodName << " Position located "
702 message << "Position located " << 710 << ( inMother == kSurface ? " on Surface " : " outside " )
703 << ( inMother == kSurface ? " on << 711 << "expected volume" << G4endl;
704 << "expected volume" << G4endl << 712 G4double safetyFromOut= motherSolid->DistanceToIn(localPosition);
705 << " Safety (from Outside) = " << 713 G4cerr << " Safety (from Outside) = " << safetyFromOut << G4endl;
706 << motherSolid->DistanceToIn(loc <<
707 G4Exception("G4VIntersectionLocator::Loc <<
708 "GeomNav1002", JustWarning, <<
709 } 714 }
710 715
711 // 1. Simple next step - quick relocation 716 // 1. Simple next step - quick relocation and check result.
712 // nav->LocateGlobalPointWithinVolume( pos 717 // nav->LocateGlobalPointWithinVolume( position );
713 718
714 // 2. Full relocation - to cross-check ans 719 // 2. Full relocation - to cross-check answer !
715 G4VPhysicalVolume* nextPhysical= nav->Loca 720 G4VPhysicalVolume* nextPhysical= nav->LocateGlobalPointAndSetup(position);
716 if( (nextPhysical != motherPhys) 721 if( (nextPhysical != motherPhys)
717 || (nextPhysical->GetCopyNo() != mothe 722 || (nextPhysical->GetCopyNo() != motherCopyNo )
718 ) 723 )
719 { 724 {
720 G4Exception("G4VIntersectionLocator::Loc << 725 G4cerr << " ERROR in " << MethodName
721 "GeomNav1002", JustWarning, << 726 << " Position located outside expected volume" << G4endl;
722 "Position located outside ex <<
723 } 727 }
724 nav->CheckMode(navCheck); // Recover orig 728 nav->CheckMode(navCheck); // Recover original value
725 } 729 }
726 else 730 else
727 { 731 {
728 nav->LocateGlobalPointWithinVolume( positi 732 nav->LocateGlobalPointWithinVolume( position );
729 } 733 }
730 return good; 734 return good;
731 } 735 }
732 736
733 ////////////////////////////////////////////// 737 ///////////////////////////////////////////////////////////////////////////
734 // 738 //
735 // LocateGlobalPointWithinVolumeCheckAndReport 739 // LocateGlobalPointWithinVolumeCheckAndReport.
736 // 740 //
737 void G4VIntersectionLocator:: 741 void G4VIntersectionLocator::
738 LocateGlobalPointWithinVolumeCheckAndReport( c 742 LocateGlobalPointWithinVolumeCheckAndReport( const G4ThreeVector& position,
739 c 743 const G4String& CodeLocationInfo,
740 G << 744 G4int CheckMode )
741 { 745 {
742 // Save value of Check mode first 746 // Save value of Check mode first
743 G4bool oldCheck = GetCheckMode(); << 747 G4bool oldCheck= GetCheckMode();
744 748
745 G4bool ok = LocateGlobalPointWithinVolumeAnd << 749 G4bool ok= LocateGlobalPointWithinVolumeAndCheck( position );
746 if( !ok ) 750 if( !ok )
747 { 751 {
748 std::ostringstream message; << 752 G4cerr << " ERROR occured in Intersection Locator" << G4endl;
749 message << "Failed point location." << G4e << 753 G4cerr << " Code Location info: " << CodeLocationInfo << G4endl;
750 << " Code Location info: " << Co << 754 if( CheckMode > 1 ) {
751 G4Exception("G4VIntersectionLocator::Locat << 755 // Additional information
752 "GeomNav1002", JustWarning, me << 756
>> 757 }
753 } 758 }
754 759
755 SetCheckMode( oldCheck ); 760 SetCheckMode( oldCheck );
756 } 761 }
757 762
758 ////////////////////////////////////////////// 763 ///////////////////////////////////////////////////////////////////////////
759 // 764 //
760 // ReportReversedPoints. 765 // ReportReversedPoints.
761 // 766 //
762 void G4VIntersectionLocator:: 767 void G4VIntersectionLocator::
763 ReportReversedPoints( std::ostringstream& msg, 768 ReportReversedPoints( std::ostringstream& msg,
764 const G4FieldTrack& Star 769 const G4FieldTrack& StartPointVel,
765 const G4FieldTrack& EndP 770 const G4FieldTrack& EndPointVel,
766 G4double NewSafety 771 G4double NewSafety, G4double epsStep,
767 const G4FieldTrack& A_Pt 772 const G4FieldTrack& A_PtVel,
768 const G4FieldTrack& B_Pt 773 const G4FieldTrack& B_PtVel,
769 const G4FieldTrack& SubS 774 const G4FieldTrack& SubStart_PtVel,
770 const G4ThreeVector& E_P 775 const G4ThreeVector& E_Point,
771 const G4FieldTrack& Appr 776 const G4FieldTrack& ApproxIntersecPointV,
772 G4int substep_no, 777 G4int substep_no, G4int substep_no_p, G4int depth )
773 { 778 {
774 // Expect that 'msg' can hold the name of t 779 // Expect that 'msg' can hold the name of the calling method
775 780
776 // FieldTrack 'points' A and B have been ta 781 // FieldTrack 'points' A and B have been tangled
777 // Whereas A should be before B, it is foun 782 // Whereas A should be before B, it is found that curveLen(B) < curveLen(A)
778 G4int verboseLevel= 5; 783 G4int verboseLevel= 5;
779 G4double curveDist = B_PtVel.GetCurveLength 784 G4double curveDist = B_PtVel.GetCurveLength() - A_PtVel.GetCurveLength();
780 G4VIntersectionLocator::printStatus( A_PtVe 785 G4VIntersectionLocator::printStatus( A_PtVel, B_PtVel,
781 -1.0, NewSafety, s 786 -1.0, NewSafety, substep_no, msg, verboseLevel );
782 msg << "Error in advancing propagation." << 787 msg << "Error in advancing propagation." << G4endl
783 << " The final curve point is NOT fur << 788 << " Point A (start) is " << A_PtVel << G4endl
784 << " than the original!" << G4endl << 789 << " Point B (end) is " << B_PtVel << G4endl
785 << " Going *backwards* from len(A) = << 790 << " Curve distance is " << curveDist << G4endl
786 << " to len(B) = " << B_PtVel.GetCurve <<
787 << " Curve distance is " << curveD <<
788 << G4endl 791 << G4endl
789 << " Point A' (start) is " << A_Pt << 792 << "The final curve point is not further along"
790 << " Point B' (end) is " << B_Pt << 793 << " than the original!" << G4endl;
791 msg << " fEpsStep= " << epsStep << G4e << 794 msg << " Value of fEpsStep= " << epsStep << G4endl;
792 << 795
793 G4long oldprc = msg.precision(20); << 796 G4int oldprc = msg.precision(20);
794 msg << " In full precision, the position, m << 797 msg << " Point A (Curve start) is " << StartPointVel << G4endl
795 << " ... are: " << G4endl; << 798 << " Point B (Curve end) is " << EndPointVel << G4endl
796 msg << " Point A[0] (Curve start) is " << << 799 << " Point A (Current start) is " << A_PtVel << G4endl
797 << " Point S (Sub start) is " << << 800 << " Point B (Current end) is " << B_PtVel << G4endl
798 << " Point A' (Current start) is " << << 801 << " Point S (Sub start) is " << SubStart_PtVel
799 << " Point E (Trial Point) is " << << 802 << " Point E (Trial Point) is " << E_Point << G4endl
800 << " Point F (Intersection) is " << << 803 << " Point F (Intersection) is " << ApproxIntersecPointV
801 << " Point B' (Current end) is " << <<
802 << " Point B[0] (Curve end) is " << <<
803 << G4endl 804 << G4endl
804 << " LocateIntersection parameters are 805 << " LocateIntersection parameters are : " << G4endl
805 << " Substep no (total) = " << su 806 << " Substep no (total) = " << substep_no << G4endl
806 << " Substep no = " << su << 807 << " Substep (depth= " << depth << substep_no_p;
807 msg.precision(oldprc); 808 msg.precision(oldprc);
808 } 809 }
809 810
810 ////////////////////////////////////////////// 811 ///////////////////////////////////////////////////////////////////////////
811 // 812 //
812 // ReportProgress. 813 // ReportProgress.
813 // 814 //
814 void G4VIntersectionLocator::ReportProgress( s 815 void G4VIntersectionLocator::ReportProgress( std::ostream& oss,
815 const G4Fi 816 const G4FieldTrack& StartPointVel,
816 const G4Fi 817 const G4FieldTrack& EndPointVel,
817 G4in 818 G4int substep_no,
818 const G4Fi 819 const G4FieldTrack& A_PtVel,
819 const G4Fi 820 const G4FieldTrack& B_PtVel,
820 G4do 821 G4double safetyLast,
821 G4in 822 G4int depth )
822 823
823 { 824 {
824 oss << "ReportProgress: Current status of in 825 oss << "ReportProgress: Current status of intersection search: " << G4endl;
825 if( depth > 0 ) { oss << " Depth= " << depth << 826 if( depth > 0 ) oss << " Depth= " << depth;
826 oss << " Substep no = " << substep_no << G4e 827 oss << " Substep no = " << substep_no << G4endl;
827 G4int verboseLevel = 5; 828 G4int verboseLevel = 5;
>> 829
>> 830 // printStatus args: (FT0, FT1, dRequestStep, dSafety, iStepNum, os, iVerb);
828 G4double safetyPrev = -1.0; // Add as argum 831 G4double safetyPrev = -1.0; // Add as argument ?
829 832
830 printStatus( StartPointVel, EndPointVel, -1. 833 printStatus( StartPointVel, EndPointVel, -1.0, -1.0, -1,
831 oss, verboseLevel); 834 oss, verboseLevel);
832 oss << " * Start and end-point of requested 835 oss << " * Start and end-point of requested Step:" << G4endl;
833 oss << " ** State of point A: "; 836 oss << " ** State of point A: ";
834 printStatus( A_PtVel, A_PtVel, -1.0, safetyP 837 printStatus( A_PtVel, A_PtVel, -1.0, safetyPrev, substep_no-1,
835 oss, verboseLevel); 838 oss, verboseLevel);
836 oss << " ** State of point B: "; 839 oss << " ** State of point B: ";
837 printStatus( A_PtVel, B_PtVel, -1.0, safetyL 840 printStatus( A_PtVel, B_PtVel, -1.0, safetyLast, substep_no,
838 oss, verboseLevel); 841 oss, verboseLevel);
839 } 842 }
840 843
841 ////////////////////////////////////////////// 844 ///////////////////////////////////////////////////////////////////////////
842 // 845 //
843 // ReportImmediateHit. 846 // ReportImmediateHit.
844 // 847 //
845 void 848 void
846 G4VIntersectionLocator::ReportImmediateHit( co 849 G4VIntersectionLocator::ReportImmediateHit( const char* MethodName,
847 co 850 const G4ThreeVector& StartPosition,
848 co 851 const G4ThreeVector& TrialPoint,
849 852 G4double tolerance,
850 un 853 unsigned long int numCalls )
851 { 854 {
852 static G4ThreadLocal unsigned int occurredOn 855 static G4ThreadLocal unsigned int occurredOnTop= 0;
853 static G4ThreadLocal G4ThreeVector* ptrLast << 856 static G4ThreadLocal G4ThreeVector *ptrLast= 0;
854 if( ptrLast == nullptr ) << 857 if( !ptrLast )
855 { 858 {
856 ptrLast= new G4ThreeVector( DBL_MAX, DBL_ 859 ptrLast= new G4ThreeVector( DBL_MAX, DBL_MAX, DBL_MAX );
857 G4AutoDelete::Register(ptrLast); 860 G4AutoDelete::Register(ptrLast);
858 } 861 }
859 G4ThreeVector &lastStart= *ptrLast; 862 G4ThreeVector &lastStart= *ptrLast;
860 863
861 if( (TrialPoint - StartPosition).mag2() < to 864 if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
862 { 865 {
863 static G4ThreadLocal unsigned int numUnmo << 866 static G4ThreadLocal unsigned int numUnmoved= 0;
864 static G4ThreadLocal unsigned int numStil << 867 static G4ThreadLocal unsigned int numStill= 0; // Still at same point
865 868
866 G4cout << "Intersection F == start A in " 869 G4cout << "Intersection F == start A in " << MethodName;
867 G4cout << "Start Point: " << StartPositio 870 G4cout << "Start Point: " << StartPosition << G4endl;
868 G4cout << " Start-Trial: " << TrialPoint << 871 // G4cout << "Trial Point: " << TrialPoint << G4endl;
>> 872 G4cout << " Start-Trial: " << TrialPoint - StartPosition; // << G4endl;
869 G4cout << " Start-last: " << StartPositio 873 G4cout << " Start-last: " << StartPosition - lastStart;
870 874
871 if( (StartPosition - lastStart).mag() < t 875 if( (StartPosition - lastStart).mag() < tolerance )
872 { 876 {
873 // We are at position of last 'Start' 877 // We are at position of last 'Start' position - ie unmoved
874 ++numUnmoved; 878 ++numUnmoved;
875 ++numStill; 879 ++numStill;
876 G4cout << " { Unmoved: " << " still#= 880 G4cout << " { Unmoved: " << " still#= " << numStill
877 << " total # = " << numUnmoved 881 << " total # = " << numUnmoved << " } - ";
878 } 882 }
879 else 883 else
880 { 884 {
881 numStill = 0; 885 numStill = 0;
882 } 886 }
883 G4cout << " Occurred: " << ++occurredOnTo << 887 G4cout << " Occured: " << ++occurredOnTop;
884 G4cout << " out of total calls= " << num 888 G4cout << " out of total calls= " << numCalls;
885 G4cout << G4endl; 889 G4cout << G4endl;
886 lastStart = StartPosition; 890 lastStart = StartPosition;
887 } 891 }
888 } // End of ReportImmediateHit() 892 } // End of ReportImmediateHit()
889 893