source: trunk/source/geometry/navigation/src/G4VIntersectionLocator.cc@ 1317

Last change on this file since 1317 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 13.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4VIntersectionLocator.cc,v 1.7 2009/11/27 15:21:59 japost Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29// Class G4VIntersectionLocator implementation
30//
31// 27.10.08 - John Apostolakis, Tatiana Nikitina.
32// ---------------------------------------------------------------------------
33
34#include <iomanip>
35
36#include "globals.hh"
37#include "G4ios.hh"
38#include "G4VIntersectionLocator.hh"
39#include "G4GeometryTolerance.hh"
40
41///////////////////////////////////////////////////////////////////////////
42//
43// Constructor
44//
45G4VIntersectionLocator:: G4VIntersectionLocator(G4Navigator *theNavigator):
46 fUseNormalCorrection(false),
47 fiNavigator( theNavigator ),
48 fiChordFinder( 0 ), // Not set - overridden at each step
49 fiEpsilonStep( -1.0 ), // Out of range - overridden at each step
50 fiDeltaIntersection( -1.0 ), // Out of range - overridden at each step
51 fiUseSafety(false) // Default - overridden at each step
52{
53 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
54 fVerboseLevel = 0;
55 fHelpingNavigator = new G4Navigator();
56}
57
58///////////////////////////////////////////////////////////////////////////
59//
60// Destructor.
61//
62G4VIntersectionLocator::~G4VIntersectionLocator()
63{
64 delete fHelpingNavigator;
65}
66
67///////////////////////////////////////////////////////////////////////////
68//
69// Dumps status of propagator.
70//
71void
72G4VIntersectionLocator::printStatus( const G4FieldTrack& StartFT,
73 const G4FieldTrack& CurrentFT,
74 G4double requestStep,
75 G4double safety,
76 G4int stepNo)
77{
78 const G4int verboseLevel= fVerboseLevel;
79 const G4ThreeVector StartPosition = StartFT.GetPosition();
80 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
81 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
82 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
83
84 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
85
86 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
87 {
88 static G4int noPrecision= 4;
89 G4cout.precision(noPrecision);
90 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
91 G4cout << std::setw( 6) << " "
92 << std::setw( 25) << " Current Position and Direction" << " "
93 << G4endl;
94 G4cout << std::setw( 5) << "Step#"
95 << std::setw(10) << " s " << " "
96 << std::setw(10) << "X(mm)" << " "
97 << std::setw(10) << "Y(mm)" << " "
98 << std::setw(10) << "Z(mm)" << " "
99 << std::setw( 7) << " N_x " << " "
100 << std::setw( 7) << " N_y " << " "
101 << std::setw( 7) << " N_z " << " " ;
102 // << G4endl;
103 G4cout // << " >>> "
104 << std::setw( 7) << " Delta|N|" << " "
105 // << std::setw( 7) << " Delta(N_z) " << " "
106 << std::setw( 9) << "StepLen" << " "
107 << std::setw(12) << "StartSafety" << " "
108 << std::setw( 9) << "PhsStep" << " ";
109
110 G4cout << G4endl;
111 }
112 if((stepNo == 0) && (verboseLevel <=3))
113 {
114 // Recurse to print the start values
115 //
116 printStatus( StartFT, StartFT, -1.0, safety, -1);
117 }
118 if( verboseLevel <= 3 )
119 {
120 if( stepNo >= 0)
121 {
122 G4cout << std::setw( 4) << stepNo << " ";
123 }
124 else
125 {
126 G4cout << std::setw( 5) << "Start" ;
127 }
128 G4cout.precision(8);
129 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
130 G4cout.precision(8);
131 G4cout << std::setw(10) << CurrentPosition.x() << " "
132 << std::setw(10) << CurrentPosition.y() << " "
133 << std::setw(10) << CurrentPosition.z() << " ";
134 G4cout.precision(4);
135 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
136 << std::setw( 7) << CurrentUnitVelocity.y() << " "
137 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
138 // G4cout << G4endl;
139 // G4cout << " >>> " ;
140 G4cout.precision(3);
141 G4cout << std::setw( 7)
142 << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag()
143 << " ";
144 // << std::setw( 7)
145 // << CurrentUnitVelocity.z() - InitialUnitVelocity.z() << " ";
146 G4cout << std::setw( 9) << step_len << " ";
147 G4cout << std::setw(12) << safety << " ";
148 if( requestStep != -1.0 )
149 {
150 G4cout << std::setw( 9) << requestStep << " ";
151 }
152 else
153 {
154 G4cout << std::setw( 9) << "Init/NotKnown" << " ";
155 }
156 G4cout << G4endl;
157 }
158 else // if( verboseLevel > 3 )
159 {
160 // Multi-line output
161
162 G4cout << "Step taken was " << step_len
163 << " out of PhysicalStep= " << requestStep << G4endl;
164 G4cout << "Final safety is: " << safety << G4endl;
165
166 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
167 << G4endl;
168 G4cout << G4endl;
169 }
170}
171
172///////////////////////////////////////////////////////////////////////////
173//
174// ReEstimateEndPoint.
175//
176G4FieldTrack G4VIntersectionLocator::
177ReEstimateEndpoint( const G4FieldTrack &CurrentStateA,
178 const G4FieldTrack &EstimatedEndStateB,
179 G4double linearDistSq,
180 G4double curveDist )
181{
182 G4FieldTrack newEndPoint( CurrentStateA );
183 G4MagInt_Driver* integrDriver= GetChordFinderFor()->GetIntegrationDriver();
184
185 G4FieldTrack retEndPoint( CurrentStateA );
186 G4bool goodAdvance;
187 G4int itrial=0;
188 const G4int no_trials= 20;
189
190 G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
191 do
192 {
193 G4double currentCurveLen= newEndPoint.GetCurveLength();
194 G4double advanceLength= endCurveLen - currentCurveLen ;
195 if (std::abs(advanceLength)<kCarTolerance)
196 {
197 goodAdvance=true;
198 }
199 else{
200 goodAdvance=
201 integrDriver->AccurateAdvance(newEndPoint, advanceLength,
202 GetEpsilonStepFor());
203 }
204 }
205 while( !goodAdvance && (++itrial < no_trials) );
206
207 if( goodAdvance )
208 {
209 retEndPoint= newEndPoint;
210 }
211 else
212 {
213 retEndPoint= EstimatedEndStateB; // Could not improve without major work !!
214 }
215
216 // All the work is done
217 // below are some diagnostics only -- before the return!
218 //
219 static const G4String MethodName("G4VIntersectionLocator::ReEstimateEndpoint");
220
221#ifdef G4VERBOSE
222 G4int latest_good_trials=0;
223 if( itrial > 1)
224 {
225 if( fVerboseLevel > 0 )
226 {
227 G4cout << MethodName << " called - goodAdv= " << goodAdvance
228 << " trials = " << itrial
229 << " previous good= " << latest_good_trials
230 << G4endl;
231 }
232 latest_good_trials=0;
233 }
234 else
235 {
236 latest_good_trials++;
237 }
238#endif
239
240#ifdef G4DEBUG_FIELD
241 G4double lengthDone = newEndPoint.GetCurveLength()
242 - CurrentStateA.GetCurveLength();
243 if( !goodAdvance )
244 {
245 if( fVerboseLevel >= 3 )
246 {
247 G4cout << MethodName << "> AccurateAdvance failed " ;
248 G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
249 G4cout << " It went only " << lengthDone << " instead of " << curveDist
250 << " -- a difference of " << curveDist - lengthDone << G4endl;
251 G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
252 << G4endl;
253 }
254 }
255
256 static G4int noInaccuracyWarnings = 0;
257 G4int maxNoWarnings = 10;
258 if ( (noInaccuracyWarnings < maxNoWarnings )
259 || (fVerboseLevel > 1) )
260 {
261 G4cerr << "G4PropagatorInField::LocateIntersectionPoint():"
262 << G4endl
263 << " Warning: Integration inaccuracy requires"
264 << " an adjustment in the step's endpoint." << G4endl
265 << " Two mid-points are further apart than their"
266 << " curve length difference" << G4endl
267 << " Dist = " << std::sqrt(linearDistSq)
268 << " curve length = " << curveDist << G4endl;
269 G4cerr << " Correction applied is "
270 << (newEndPoint.GetPosition()-EstimatedEndStateB.GetPosition()).mag()
271 << G4endl;
272 }
273#else
274 // Statistics on the RMS value of the corrections
275
276 static G4int noCorrections=0;
277 static G4double sumCorrectionsSq = 0;
278 noCorrections++;
279 if( goodAdvance )
280 {
281 sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
282 newEndPoint.GetPosition()).mag2();
283 }
284 linearDistSq -= curveDist; // To use linearDistSq ... !
285#endif
286
287 return retEndPoint;
288}
289
290///////////////////////////////////////////////////////////////////////////
291//
292// Method for finding SurfaceNormal of Intersecting Solid
293//
294G4ThreeVector G4VIntersectionLocator::
295GetLocalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
296{
297 G4ThreeVector Normal(G4ThreeVector(0,0,0));
298 G4VPhysicalVolume* located;
299
300 validNormal = false;
301 fHelpingNavigator->SetWorldVolume(GetNavigatorFor()->GetWorldVolume());
302 located = fHelpingNavigator->LocateGlobalPointAndSetup( CurrentE_Point );
303 G4TouchableHistoryHandle aTouchable = fHelpingNavigator
304 ->CreateTouchableHistoryHandle();
305 G4ThreeVector localPosition = aTouchable->GetHistory()
306 ->GetTopTransform().TransformPoint(CurrentE_Point);
307
308 if( located != 0)
309 {
310 if (located->GetLogicalVolume()
311 ->GetSolid()->Inside(localPosition)==kSurface)
312 {
313 Normal = located->GetLogicalVolume()
314 ->GetSolid()->SurfaceNormal(localPosition);
315 validNormal = true;
316 }
317 }
318
319 return Normal;
320}
321
322///////////////////////////////////////////////////////////////////////////
323//
324// Adjustment of Found Intersection
325//
326G4bool G4VIntersectionLocator::
327AdjustmentOfFoundIntersection( const G4ThreeVector &CurrentA_Point,
328 const G4ThreeVector &CurrentE_Point,
329 const G4ThreeVector &CurrentF_Point,
330 const G4ThreeVector &MomentumDir,
331 const G4bool IntersectAF,
332 G4ThreeVector &IntersectionPoint, // I/O
333 G4double &NewSafety, // I/O
334 G4double &fPreviousSafety, // I/O
335 G4ThreeVector &fPreviousSftOrigin )// I/O
336{
337 G4double dist,lambda;
338 G4ThreeVector Normal, NewPoint, Point_G;
339 G4bool goodAdjust=false, Intersects_FP=false, validNormal=false;
340
341 // Get SurfaceNormal of Intersecting Solid
342 //
343 Normal=GetLocalSurfaceNormal(CurrentE_Point,validNormal);
344 if(!validNormal) { return false; }
345
346 // Intersection between Line and Plane
347 //
348 G4double n_d_m = Normal.dot(MomentumDir);
349 if ( std::abs(n_d_m)>kCarTolerance )
350 {
351 if ( fVerboseLevel>1 )
352 {
353 G4cerr << "WARNING - "
354 << "G4VIntersectionLocator::AdjustementOfFoundIntersection()"
355 << G4endl
356 << " No intersection. Parallels lines!" << G4endl;
357 return false;
358 }
359 lambda =- Normal.dot(CurrentF_Point-CurrentE_Point)/n_d_m;
360
361 // New candidate for Intersection
362 //
363 NewPoint = CurrentF_Point+lambda*MomentumDir;
364
365 // Distance from CurrentF to Calculated Intersection
366 //
367 dist = std::abs(lambda);
368
369 if ( dist<kCarTolerance*0.001 ) { return false; }
370
371 // Calculation of new intersection point on the path.
372 //
373 if ( IntersectAF ) // First part intersects
374 {
375 G4double stepLengthFP;
376 G4ThreeVector Point_P = CurrentA_Point;
377 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_P);
378 Intersects_FP = IntersectChord( Point_P, NewPoint, NewSafety,
379 fPreviousSafety, fPreviousSftOrigin,
380 stepLengthFP, Point_G );
381
382 }
383 else // Second part intersects
384 {
385 G4double stepLengthFP;
386 GetNavigatorFor()->LocateGlobalPointWithinVolume(CurrentF_Point );
387 Intersects_FP = IntersectChord( CurrentF_Point, NewPoint, NewSafety,
388 fPreviousSafety, fPreviousSftOrigin,
389 stepLengthFP, Point_G );
390 }
391 if ( Intersects_FP )
392 {
393 goodAdjust = true;
394 IntersectionPoint = Point_G;
395 }
396 }
397
398 return goodAdjust;
399}
Note: See TracBrowser for help on using the repository browser.