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 | // |
---|
45 | G4VIntersectionLocator:: 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 | // |
---|
62 | G4VIntersectionLocator::~G4VIntersectionLocator() |
---|
63 | { |
---|
64 | delete fHelpingNavigator; |
---|
65 | } |
---|
66 | |
---|
67 | /////////////////////////////////////////////////////////////////////////// |
---|
68 | // |
---|
69 | // Dumps status of propagator. |
---|
70 | // |
---|
71 | void |
---|
72 | G4VIntersectionLocator::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 | // |
---|
176 | G4FieldTrack G4VIntersectionLocator:: |
---|
177 | ReEstimateEndpoint( 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 | // |
---|
294 | G4ThreeVector G4VIntersectionLocator:: |
---|
295 | GetLocalSurfaceNormal(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 | // |
---|
326 | G4bool G4VIntersectionLocator:: |
---|
327 | AdjustmentOfFoundIntersection( 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 | } |
---|