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

Last change on this file since 1158 was 1058, checked in by garnier, 17 years ago

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