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