source: trunk/source/geometry/navigation/src/G4PropagatorInField.cc@ 857

Last change on this file since 857 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 57.0 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//
27// $Id: G4PropagatorInField.cc,v 1.43 2008/05/28 09:12:23 tnikitin Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31// This class implements an algorithm to track a particle in a
32// non-uniform magnetic field. It utilises an ODE solver (with
33// the Runge - Kutta method) to evolve the particle, and drives it
34// until the particle has traveled a set distance or it enters a new
35// volume.
36//
37// 14.10.96 John Apostolakis, design and implementation
38// 17.03.97 John Apostolakis, renaming new set functions being added
39//
40// ---------------------------------------------------------------------------
41
42#include "G4PropagatorInField.hh"
43#include "G4ios.hh"
44#include <iomanip>
45
46#include "G4ThreeVector.hh"
47#include "G4VPhysicalVolume.hh"
48#include "G4Navigator.hh"
49#include "G4GeometryTolerance.hh"
50#include "G4VCurvedTrajectoryFilter.hh"
51#include "G4ChordFinder.hh"
52
53///////////////////////////////////////////////////////////////////////////
54//
55// Constructors and destructor
56
57G4PropagatorInField::G4PropagatorInField( G4Navigator *theNavigator,
58 G4FieldManager *detectorFieldMgr )
59 : fDetectorFieldMgr(detectorFieldMgr),
60 fCurrentFieldMgr(detectorFieldMgr),
61 fNavigator(theNavigator),
62 End_PointAndTangent(G4ThreeVector(0.,0.,0.),
63 G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0),
64 fParticleIsLooping(false),
65 fVerboseLevel(0),
66 fMax_loop_count(1000),
67 fNoZeroStep(0),
68 fCharge(0.0), fInitialMomentumModulus(0.0), fMass(0.0),
69 fUseSafetyForOptimisation(true), // (false) is less sensitive to incorrect safety
70 fSetFieldMgr(false),
71 fpTrajectoryFilter( 0 )
72{
73 if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
74 else { fEpsilonStep= 1.0e-5; }
75 fActionThreshold_NoZeroSteps = 2;
76 fSevereActionThreshold_NoZeroSteps = 10;
77 fAbandonThreshold_NoZeroSteps = 50;
78 fFull_CurveLen_of_LastAttempt = -1;
79 fLast_ProposedStepLength = -1;
80 fLargestAcceptableStep = 1000.0 * meter;
81
82 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
83 fPreviousSafety= 0.0;
84 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
85 //
86 fUseBrentLocator=true;
87 // In case of too slow progress in finding Intersection Point
88 // intermediates Points on the Track must be stored.
89 // Initialise the array of Pointers [max_depth+1] to do this
90 G4ThreeVector zeroV(0.0,0.0,0.0);
91 for (G4int idepth=0; idepth<max_depth+1; idepth++ )
92 {
93 ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
94 }
95 // Counter for Maximum Number Of Trial before Intersection Found
96 maxNumberOfStepsForIntersection=0;
97 // Counter for Number Of Calls to ReIntegrationEndPoint Method
98 maxNumberOfCallsToReIntegration=0;
99 maxNumberOfCallsToReIntegration_depth=0;
100}
101
102G4PropagatorInField::~G4PropagatorInField()
103{
104 for ( G4int idepth=0; idepth<max_depth+1; idepth++)
105 {
106 delete ptrInterMedFT[idepth];
107 }
108 if(fVerboseLevel>0){
109 G4cout<<"G4PropagatorInField::Location with Max Number of Steps="
110 << maxNumberOfStepsForIntersection<<G4endl;
111 G4cout<<"G4PropagatorInField::ReIntegrateEndPoint was called "<<
112 maxNumberOfCallsToReIntegration<<" times and for depth algorithm "<<
113 maxNumberOfCallsToReIntegration_depth<<" times"<<G4endl;
114 }
115}
116
117///////////////////////////////////////////////////////////////////////////
118//
119// Compute the next geometric Step
120
121G4double
122G4PropagatorInField::ComputeStep(
123 G4FieldTrack& pFieldTrack,
124 G4double CurrentProposedStepLength,
125 G4double& currentSafety, // IN/OUT
126 G4VPhysicalVolume* pPhysVol)
127{
128
129 // If CurrentProposedStepLength is too small for finding Chords
130 // then return with no action (for now - TODO: some action)
131 //
132 if(CurrentProposedStepLength<kCarTolerance)
133 {
134 return kInfinity;
135 }
136
137 // Introducing smooth trajectory display (jacek 01/11/2002)
138 //
139 if (fpTrajectoryFilter)
140 {
141 fpTrajectoryFilter->CreateNewTrajectorySegment();
142 }
143
144 // Parameters for adaptive Runge-Kutta integration
145
146 G4double h_TrialStepSize; // 1st Step Size
147 G4double TruePathLength = CurrentProposedStepLength;
148 G4double StepTaken = 0.0;
149 G4double s_length_taken, epsilon ;
150 G4bool intersects;
151 G4bool first_substep = true;
152
153 G4double NewSafety;
154 fParticleIsLooping = false;
155
156 // If not yet done,
157 // Set the field manager to the local one if the volume has one,
158 // or to the global one if not
159 //
160 if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol );
161 // For the next call, the field manager must again be set
162 fSetFieldMgr= false;
163
164 GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass);
165
166 G4FieldTrack CurrentState(pFieldTrack);
167 G4FieldTrack OriginalState = CurrentState;
168
169 // If the Step length is "infinite", then an approximate-maximum Step
170 // length (used to calculate the relative accuracy) must be guessed.
171 //
172 if( CurrentProposedStepLength >= fLargestAcceptableStep )
173 {
174 G4ThreeVector StartPointA, VelocityUnit;
175 StartPointA = pFieldTrack.GetPosition();
176 VelocityUnit = pFieldTrack.GetMomentumDir();
177
178 G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
179 fNavigator->GetWorldVolume()->GetLogicalVolume()->
180 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
181 CurrentProposedStepLength= std::min( trialProposedStep,
182 fLargestAcceptableStep );
183 }
184 epsilon = GetDeltaOneStep() / CurrentProposedStepLength;
185 // G4double raw_epsilon= epsilon;
186 G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
187 G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();;
188 if( epsilon < epsilonMin ) epsilon = epsilonMin;
189 if( epsilon > epsilonMax ) epsilon = epsilonMax;
190 SetEpsilonStep( epsilon );
191
192 // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
193 // << " final= " << epsilon << G4endl;
194
195 // Shorten the proposed step in case of earlier problems (zero steps)
196 //
197 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
198 {
199 G4double stepTrial;
200
201 stepTrial= fFull_CurveLen_of_LastAttempt;
202 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
203 stepTrial= fLast_ProposedStepLength;
204
205 G4double decreaseFactor = 0.9; // Unused default
206 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
207 && (stepTrial > 1000.0*kCarTolerance) )
208 {
209 // Ensure quicker convergence
210 //
211 decreaseFactor= 0.1;
212 }
213 else
214 {
215 // We are in significant difficulties, probably at a boundary that
216 // is either geometrically sharp or between very different materials.
217 // Careful decreases to cope with tolerance are required.
218 //
219 if( stepTrial > 1000.0*kCarTolerance )
220 decreaseFactor = 0.25; // Try slow decreases
221 else if( stepTrial > 100.0*kCarTolerance )
222 decreaseFactor= 0.5; // Try slower decreases
223 else if( stepTrial > 10.0*kCarTolerance )
224 decreaseFactor= 0.75; // Try even slower decreases
225 else
226 decreaseFactor= 0.9; // Try very slow decreases
227 }
228 stepTrial *= decreaseFactor;
229
230#ifdef G4DEBUG_FIELD
231 PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
232 stepTrial, pFieldTrack);
233#endif
234 if( stepTrial == 0.0 )
235 {
236 G4cout << " G4PropagatorInField::ComputeStep "
237 << " Particle abandoned due to lack of progress in field."
238 << G4endl
239 << " Properties : " << pFieldTrack << " "
240 << G4endl;
241 G4cerr << " G4PropagatorInField::ComputeStep "
242 << " ERROR : attempting a zero step= " << stepTrial << G4endl
243 << " while attempting to progress after " << fNoZeroStep
244 << " trial steps. Will abandon step." << G4endl;
245 fParticleIsLooping= true;
246 return 0; // = stepTrial;
247 }
248 if( stepTrial < CurrentProposedStepLength )
249 CurrentProposedStepLength = stepTrial;
250 }
251 fLast_ProposedStepLength = CurrentProposedStepLength;
252
253 G4int do_loop_count = 0;
254 do
255 {
256 G4FieldTrack SubStepStartState = CurrentState;
257 G4ThreeVector SubStartPoint = CurrentState.GetPosition();
258
259 if( !first_substep) {
260 fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
261 }
262
263 // How far to attempt to move the particle !
264 //
265 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
266
267 // Integrate as far as "chord miss" rule allows.
268 //
269 s_length_taken = GetChordFinder()->AdvanceChordLimited(
270 CurrentState, // Position & velocity
271 h_TrialStepSize,
272 fEpsilonStep,
273 fPreviousSftOrigin,
274 fPreviousSafety
275 );
276 // CurrentState is now updated with the final position and velocity.
277
278 fFull_CurveLen_of_LastAttempt = s_length_taken;
279
280 G4ThreeVector EndPointB = CurrentState.GetPosition();
281 G4ThreeVector InterSectionPointE;
282 G4double LinearStepLength;
283
284 // Intersect chord AB with geometry
285 intersects= IntersectChord( SubStartPoint, EndPointB,
286 NewSafety, LinearStepLength,
287 InterSectionPointE );
288 // E <- Intersection Point of chord AB and either volume A's surface
289 // or a daughter volume's surface ..
290
291 if( first_substep ) {
292 currentSafety = NewSafety;
293 } // Updating safety in other steps is potential future extention
294
295 if( intersects )
296 {
297 G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
298
299 // Find the intersection point of AB true path with the surface
300 // of vol(A), if it exists. Start with point E as first "estimate".
301 G4bool recalculatedEndPt= false;
302 G4bool found_intersection =
303 LocateIntersectionPoint( SubStepStartState, CurrentState,
304 InterSectionPointE, IntersectPointVelct_G,
305 recalculatedEndPt);
306 // G4cout<<"In Locate"<<recalculatedEndPt<<" and V"<<IntersectPointVelct_G.GetPosition()<<G4endl;
307 intersects = intersects && found_intersection;
308 if( found_intersection ) {
309 End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
310 StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
311 - OriginalState.GetCurveLength();
312 } else {
313 // intersects= false; // "Minor" chords do not intersect
314 if( recalculatedEndPt ){
315 CurrentState= IntersectPointVelct_G;
316 }
317 }
318 }
319 if( !intersects )
320 {
321 StepTaken += s_length_taken;
322 // For smooth trajectory display (jacek 01/11/2002)
323 if (fpTrajectoryFilter) {
324 fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
325 }
326 }
327 first_substep = false;
328
329#ifdef G4DEBUG_FIELD
330 if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
331 printStatus( SubStepStartState, // or OriginalState,
332 CurrentState, CurrentProposedStepLength,
333 NewSafety, do_loop_count, pPhysVol );
334 }
335#endif
336#ifdef G4VERBOSE
337 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
338 if( do_loop_count == fMax_loop_count-9 ){
339 G4cout << "G4PropagatorInField::ComputeStep "
340 << " Difficult track - taking many sub steps." << G4endl;
341 }
342 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
343 NewSafety, do_loop_count, pPhysVol );
344 }
345#endif
346
347 do_loop_count++;
348
349 } while( (!intersects )
350 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
351 && ( do_loop_count < fMax_loop_count ) );
352
353 if( do_loop_count >= fMax_loop_count )
354 {
355 fParticleIsLooping = true;
356
357 if ( fVerboseLevel > 0 ){
358 G4cout << "G4PropagateInField: Killing looping particle "
359 // << " of " << energy << " energy "
360 << " after " << do_loop_count << " field substeps "
361 << " totaling " << StepTaken / mm << " mm " ;
362 if( pPhysVol )
363 G4cout << " in the volume " << pPhysVol->GetName() ;
364 else
365 G4cout << " in unknown or null volume. " ;
366 G4cout << G4endl;
367 }
368 }
369
370 if( !intersects )
371 {
372 // Chord AB or "minor chords" do not intersect
373 // B is the endpoint Step of the current Step.
374 //
375 End_PointAndTangent = CurrentState;
376 TruePathLength = StepTaken;
377 }
378
379 // Set pFieldTrack to the return value
380 //
381 pFieldTrack = End_PointAndTangent;
382
383#ifdef G4VERBOSE
384 // Check that "s" is correct
385 //
386 if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
387 - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
388 {
389 G4cerr << " ERROR - G4PropagatorInField::ComputeStep():" << G4endl
390 << " Curve length mis-match, is advancement wrong ? " << G4endl;
391 G4cerr << " The curve length of the endpoint should be: "
392 << OriginalState.GetCurveLength() + TruePathLength << G4endl
393 << " and it is instead: "
394 << End_PointAndTangent.GetCurveLength() << "." << G4endl
395 << " A difference of: "
396 << OriginalState.GetCurveLength() + TruePathLength
397 - End_PointAndTangent.GetCurveLength() << G4endl;
398 G4cerr << " Original state= " << OriginalState << G4endl
399 << " Proposed state= " << End_PointAndTangent << G4endl;
400 G4Exception("G4PropagatorInField::ComputeStep()", "IncorrectProposedEndPoint",
401 FatalException,
402 "Curve length mis-match between original state and proposed endpoint of propagation.");
403 }
404#endif
405
406 // In particular anomalous cases, we can get repeated zero steps
407 // In order to correct this efficiently, we identify these cases
408 // and only take corrective action when they occur.
409 //
410 if( TruePathLength < 0.5*kCarTolerance )
411 fNoZeroStep++;
412 else
413 fNoZeroStep = 0;
414
415 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps ) {
416 fParticleIsLooping = true;
417 G4cout << " WARNING - G4PropagatorInField::ComputeStep():" << G4endl
418 << " Zero progress for " << fNoZeroStep << " attempted steps."
419 << G4endl;
420 G4cout << "Proposed Step is "<<CurrentProposedStepLength <<" but Step Taken is "<< fFull_CurveLen_of_LastAttempt <<G4endl;
421 G4cout << "For Particle with Charge ="<<fCharge
422 << " Momentum="<< fInitialMomentumModulus<<" Mass="<< fMass<<G4endl;
423 if( pPhysVol )
424 G4cout << " in the volume " << pPhysVol->GetName() ;
425 else
426 G4cout << " in unknown or null volume. " ;
427 G4cout << G4endl;
428 if ( fVerboseLevel > 2 )
429 G4cout << " Particle that is stuck will be killed." << G4endl;
430 fNoZeroStep = 0;
431 }
432
433 return TruePathLength;
434}
435
436// --------------------------------------------------------------------------
437// G4bool
438// G4PropagatorInField::LocateIntersectionPoint(
439// const G4FieldTrack& CurveStartPointVelocity, // A
440// const G4FieldTrack& CurveEndPointVelocity, // B
441// const G4ThreeVector& TrialPoint, // E
442// G4FieldTrack& IntersectedOrRecalculated // Output
443// G4bool& recalculated) // Out
444// --------------------------------------------------------------------------
445//
446// Function that returns the intersection of the true path with the surface
447// of the current volume (either the external one or the inner one with one
448// of the daughters
449//
450// A = Initial point
451// B = another point
452//
453// Both A and B are assumed to be on the true path.
454//
455// E is the first point of intersection of the chord AB with
456// a volume other than A (on the surface of A or of a daughter)
457//
458// Convention of Use :
459// i) If it returns "true", then IntersectionPointVelocity is set
460// to the approximate intersection point.
461// ii) If it returns "false", no intersection was found.
462// The validity of IntersectedOrRecalculated depends on 'recalculated'
463// a) if latter is false, then IntersectedOrRecalculated is invalid.
464// b) if latter is true, then IntersectedOrRecalculated is
465// the new endpoint, due to a re-integration.
466// --------------------------------------------------------------------------
467
468G4bool
469G4PropagatorInField::LocateIntersectionPoint(
470 const G4FieldTrack& CurveStartPointVelocity, // A
471 const G4FieldTrack& CurveEndPointVelocity, // B
472 const G4ThreeVector& TrialPoint, // E
473 G4FieldTrack& IntersectedOrRecalculatedFT, // Out: point found
474 G4bool& recalculatedEndPoint) // Out:
475{
476 // Find Intersection Point ( A, B, E ) of true path AB - start at E.
477
478 G4bool found_approximate_intersection = false;
479 G4bool there_is_no_intersection = false;
480
481 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
482 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
483 G4ThreeVector CurrentE_Point = TrialPoint;
484 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
485 G4double NewSafety= -0.0;
486
487 G4bool final_section= true; // Shows whether current section is last
488 // (i.e. B=full end)
489 G4bool first_section=true;
490 recalculatedEndPoint= false;
491
492 G4bool restoredFullEndpoint= false;
493
494 G4int substep_no = 0;
495
496 // Limits for substep number
497 //
498 const G4int max_substeps= 10000; // Test 120 (old value 100 )
499 const G4int warn_substeps= 1000; // 100
500
501 // Statistics for substeps
502 //
503 static G4int max_no_seen= -1;
504 static G4int trigger_substepno_print= warn_substeps - 20 ;
505
506 //--------------------------------------------------------------------------
507 // Algoritm for the case if progress in founding intersection is too slow.
508 // Process is defined too slow if after N=param_substeps advances on the
509 // path, it will be only 'fraction_done' of the total length.
510 // In this case the remaining length is divided in two half and
511 // the loop is restarted for each half.
512 // If progress is still too slow, the division in two halfs continue
513 // until 'max_depth'.
514 //--------------------------------------------------------------------------
515 G4double count_did_len=0.;
516 G4double count_all_len=0;
517 G4int param_substeps=100;//Test value for the maximum number of substeps
518 if(!fUseBrentLocator) param_substeps=10;// Reduced value for the maximum number
519
520 const G4double fraction_done=0.3;
521
522 G4bool Second_half=false; // First half or second half of divided step
523
524 // We need to know this for the 'final_section':
525 // real 'final_section' or first half 'final_section'
526 // In algorithm it is considered that the 'Second_half' is true
527 // and it becomes false only if we are in the first-half of level
528 // depthness or if we are in the first section
529
530 G4int depth=0; // Depth counts how many subdivisions of initial step made
531
532#ifdef G4DEBUG_FIELD
533 static G4double tolerance= 1.0e-8;
534 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
535 if( (TrialPoint - StartPosition).mag() < tolerance * mm )
536 {
537 G4cerr << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
538 << G4endl
539 << " Intermediate F point is on top of starting point A."
540 << G4endl;
541 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
542 "IntersectionPointIsAtStart", JustWarning,
543 "Intersection point F is exactly at start point A." );
544 }
545#endif
546
547 // Intermediates Points on the Track = Subdivided Points must be stored.
548 // Give the initial values to 'InterMedFt'
549 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
550 //
551 *ptrInterMedFT[0] = CurveEndPointVelocity;
552 for (G4int idepth=1; idepth<max_depth+1; idepth++ )
553 {
554 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
555 }
556
557 // 'SubStartPoint' is needed to calculate the length of the divided step
558 //
559 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
560
561 do
562 {
563 G4int substep_no_p = 0;
564 G4bool sub_final_section = false; // the same as final_section,
565 // but for 'sub_section'
566 do // REPEAT param
567 {
568 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
569 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
570
571 // F = a point on true AB path close to point E
572 // (the closest if possible)
573 //
574 if((!fUseBrentLocator)||(substep_no_p==0)){
575 ApproxIntersecPointV = GetChordFinder()
576 ->ApproxCurvePointV( CurrentA_PointVelocity,
577 CurrentB_PointVelocity,
578 CurrentE_Point,
579 fEpsilonStep );
580 // The above method is the key & most intuitive part ...
581 }
582#ifdef G4DEBUG_FIELD
583 if( ApproxIntersecPointV.GetCurveLength() >
584 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
585 {
586 G4cerr << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
587 << G4endl
588 << " Intermediate F point is more advanced than"
589 << " endpoint B." << G4endl;
590 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
591 "IntermediatePointConfusion", FatalException,
592 "Intermediate F point is past end B point" );
593 }
594#endif
595
596 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
597 if(substep_no> maxNumberOfStepsForIntersection)maxNumberOfStepsForIntersection=substep_no;
598 // First check whether EF is small - then F is a good approx. point
599 // Calculate the length and direction of the chord AF
600 //
601 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
602
603 if ( ChordEF_Vector.mag2() <= sqr(GetDeltaIntersection()) )
604 {
605 found_approximate_intersection = true;
606 // Create the "point" return value
607 //
608 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
609 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
610
611 // Note: in order to return a point on the boundary,
612 // we must return E. But it is F on the curve.
613 // So we must "cheat": we are using the position at point E
614 // and the velocity at point F !!!
615 //
616 // This must limit the length we can allow for displacement!
617 }
618 else // E is NOT close enough to the curve (ie point F)
619 {
620 // Check whether any volumes are encountered by the chord AF
621 // ---------------------------------------------------------
622 // First relocate to restore any Voxel etc information
623 // in the Navigator before calling ComputeStep()
624 //
625 fNavigator->LocateGlobalPointWithinVolume( Point_A );
626
627 G4ThreeVector PointG; // Candidate intersection point
628 G4double stepLengthAF;
629 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
630 NewSafety, stepLengthAF,
631 PointG );
632 if( Intersects_AF )
633 {
634 if(fUseBrentLocator){
635
636 G4FieldTrack EndPoint=ApproxIntersecPointV;
637 ApproxIntersecPointV= GetChordFinder()->ApproxCurvePointS(
638 CurrentA_PointVelocity,CurrentB_PointVelocity,
639 CurrentE_Point,CurrentF_Point,PointG,true,fEpsilonStep);
640 CurrentB_PointVelocity = EndPoint;
641 CurrentE_Point = PointG;
642 // By moving point B, must take care if current
643 // AF has no intersection to try current FB!!
644 //
645 final_section= false;
646
647 }
648 else{
649 // G is our new Candidate for the intersection point.
650 // It replaces "E" and we will repeat the test to see if
651 // it is a good enough approximate point for us.
652 // B <- F
653 // E <- G
654
655 CurrentB_PointVelocity = ApproxIntersecPointV;
656 CurrentE_Point = PointG;
657
658 // By moving point B, must take care if current
659 // AF has no intersection to try current FB!!
660 //
661 final_section= false;
662 }
663#ifdef G4VERBOSE
664 if( fVerboseLevel > 3 )
665 {
666 G4cout << "G4PiF::LI> Investigating intermediate point"
667 << " at s=" << ApproxIntersecPointV.GetCurveLength()
668 << " on way to full s="
669 << CurveEndPointVelocity.GetCurveLength() << G4endl;
670 }
671#endif
672 }
673 else // not Intersects_AF
674 {
675 // In this case:
676 // There is NO intersection of AF with a volume boundary.
677 // We must continue the search in the segment FB!
678 //
679 fNavigator->LocateGlobalPointWithinVolume( CurrentF_Point );
680
681 G4double stepLengthFB;
682 G4ThreeVector PointH;
683
684 // Check whether any volumes are encountered by the chord FB
685 // ---------------------------------------------------------
686
687 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
688 NewSafety, stepLengthFB,
689 PointH );
690 if( Intersects_FB )
691 {
692 if(fUseBrentLocator){
693 CurrentA_PointVelocity = ApproxIntersecPointV;
694 ApproxIntersecPointV= GetChordFinder()->ApproxCurvePointS(
695 CurrentA_PointVelocity,CurrentB_PointVelocity,
696 CurrentE_Point,Point_A,PointH,false,fEpsilonStep);
697 CurrentE_Point = PointH;
698 }
699 else{
700
701 // There is an intersection of FB with a volume boundary
702 // H <- First Intersection of Chord FB
703
704 // H is our new Candidate for the intersection point.
705 // It replaces "E" and we will repeat the test to see if
706 // it is a good enough approximate point for us.
707
708 // Note that F must be in volume volA (the same as A)
709 // (otherwise AF would meet a volume boundary!)
710 // A <- F
711 // E <- H
712
713 CurrentA_PointVelocity = ApproxIntersecPointV;
714 CurrentE_Point = PointH;
715 }
716 }
717 else // not Intersects_FB
718 {
719 // There is NO intersection of FB with a volume boundary
720
721 if( final_section )
722 {
723 // If B is the original endpoint, this means that whatever
724 // volume(s) intersected the original chord, none touch the
725 // smaller chords we have used.
726 // The value of 'IntersectedOrRecalculatedFT' returned is
727 // likely not valid
728
729 // Check on real final_section or SubEndSection
730 //
731 if( ((Second_half)&&(depth==0)) || (first_section) )
732 {
733 there_is_no_intersection = true; // real final_section
734 }
735 else
736 {
737 // end of subsection, not real final section
738 // exit from the and go to the depth-1 level
739
740 substep_no_p = param_substeps+2; // exit from the loop
741
742 // but 'Second_half' is still true because we need to find
743 // the 'CurrentE_point' for the next loop
744 //
745 Second_half = true;
746 sub_final_section = true;
747
748 }
749 }
750 else
751 {
752 // We must restore the original endpoint
753
754 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
755 CurrentB_PointVelocity = CurveEndPointVelocity;
756 restoredFullEndpoint = true;
757 }
758 } // Endif (Intersects_FB)
759 } // Endif (Intersects_AF)
760
761 // Ensure that the new endpoints are not further apart in space
762 // than on the curve due to different errors in the integration
763 //
764 G4double linDistSq, curveDist;
765 linDistSq = ( CurrentB_PointVelocity.GetPosition()
766 - CurrentA_PointVelocity.GetPosition() ).mag2();
767 curveDist = CurrentB_PointVelocity.GetCurveLength()
768 - CurrentA_PointVelocity.GetCurveLength();
769 if( curveDist*curveDist*(1+2*fEpsilonStep ) < linDistSq )
770 {
771 // Re-integrate to obtain a new B
772 //
773 G4FieldTrack newEndPointFT=
774 ReEstimateEndpoint( CurrentA_PointVelocity,
775 CurrentB_PointVelocity,
776 linDistSq, // to avoid recalculation
777 curveDist );
778 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
779 CurrentB_PointVelocity = newEndPointFT;
780 maxNumberOfCallsToReIntegration= maxNumberOfCallsToReIntegration+1;
781 #ifdef G4DEBUG_FIELD
782 G4cout<<"G4PIF::Call ReIntEnd1 linD="<<std::sqrt(linDistSq)<<" curve="<<curveDist<<" n="<<substep_no<<G4endl;
783 G4cout<<"G4PIF::Call ReIntEnd2 IntersectAF="<< Intersects_AF<<" final_section="<<final_section<<G4endl;
784 #endif
785 if( (final_section)&&(Second_half)&&(depth==0) ) // real final section
786 {
787 recalculatedEndPoint = true;
788 IntersectedOrRecalculatedFT = newEndPointFT;
789 // So that we can return it, if it is the endpoint!
790 }
791
792 }
793
794 if( curveDist < 0.0 )
795 {
796 G4cerr << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
797 << G4endl
798 << " Error in advancing propagation." << G4endl;
799 fVerboseLevel = 5; // Print out a maximum of information
800 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
801 -1.0, NewSafety, substep_no, 0 );
802 G4cerr << " Point A (start) is " << CurrentA_PointVelocity
803 << G4endl;
804 G4cerr << " Point B (end) is " << CurrentB_PointVelocity
805 << G4endl;
806 G4cerr << " Curve distance is " << curveDist << G4endl;
807 G4cerr << G4endl
808 << "The final curve point is not further along"
809 << " than the original!" << G4endl;
810 if( recalculatedEndPoint )
811 {
812 G4cerr << "Recalculation of EndPoint was called with fEpsStep= "
813 << fEpsilonStep << G4endl;
814 }
815 G4cerr.precision(20);
816 G4cerr << " Point A (Curve start) is " << CurveStartPointVelocity
817 << G4endl;
818 G4cerr << " Point B (Curve end) is " << CurveEndPointVelocity
819 << G4endl;
820 G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity
821 << G4endl;
822 G4cerr << " Point B (Current end) is " << CurrentB_PointVelocity
823 << G4endl;
824 G4cerr << " Point S (Sub start) is " << SubStart_PointVelocity
825 << G4endl;
826 G4cerr << " Point E (Trial Point) is " << CurrentE_Point
827 << G4endl;
828 G4cerr << " Point F (Intersection) is " << ApproxIntersecPointV
829 << G4endl;
830 G4cerr << " LocateIntersection parameters are : Substep no= "
831 << substep_no << G4endl;
832 G4cerr << " Substep depth no= "<< substep_no_p << " Depth= "
833 << depth << G4endl;
834 G4cerr << " did_len= "<< count_did_len << " all_len= "
835 << count_all_len << G4endl;
836 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
837 "FatalError", FatalException,
838 "Error in advancing propagation.");
839 }
840
841 if(restoredFullEndpoint)
842 {
843 final_section = restoredFullEndpoint;
844 restoredFullEndpoint = false;
845 }
846 } // EndIf ( E is close enough to the curve, ie point F. )
847 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
848
849#ifdef G4DEBUG_LOCATE_INTERSECTION
850 if( substep_no >= trigger_substepno_print )
851 {
852 G4cout << "Difficulty in converging in "
853 << "G4PropagatorInField::LocateIntersectionPoint():"
854 << G4endl
855 << " Substep no = " << substep_no << G4endl;
856 if( substep_no == trigger_substepno_print )
857 {
858 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
859 -1.0, NewSafety, 0, 0);
860 }
861 G4cout << " State of point A: ";
862 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
863 -1.0, NewSafety, substep_no-1, 0);
864 G4cout << " State of point B: ";
865 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
866 -1.0, NewSafety, substep_no, 0);
867 }
868#endif
869
870 substep_no++;
871 substep_no_p++;
872
873 } while ( ( ! found_approximate_intersection )
874 && ( ! there_is_no_intersection )
875 && ( substep_no_p <= param_substeps) ); // UNTIL found or
876 // failed param substep
877 first_section = false;
878
879 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
880 {
881 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
882 - SubStart_PointVelocity.GetCurveLength());
883 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
884 - SubStart_PointVelocity.GetCurveLength());
885 count_did_len=did_len;
886 count_all_len=all_len;
887 G4double stepLengthAB;
888 G4ThreeVector PointGe;
889
890 // Check if progress is too slow and if it possible to go deeper,
891 // then halve the step if so
892 //
893 if( ( ( did_len )<fraction_done*all_len)
894 && (depth<max_depth) && (!sub_final_section) )
895 {
896
897 Second_half=false;
898 depth++;
899
900 G4double Sub_len = (all_len-did_len)/(2.);
901 G4FieldTrack start = CurrentA_PointVelocity;
902 G4MagInt_Driver* integrDriver=GetChordFinder()->GetIntegrationDriver();
903 integrDriver->AccurateAdvance(start, Sub_len, fEpsilonStep);
904 *ptrInterMedFT[depth] = start;
905 CurrentB_PointVelocity = *ptrInterMedFT[depth];
906
907 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
908 //
909 SubStart_PointVelocity = CurrentA_PointVelocity;
910
911 // Find new trial intersection point needed at start of the loop
912 //
913 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
914 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
915
916 fNavigator->LocateGlobalPointWithinVolume(Point_A);
917 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
918 NewSafety, stepLengthAB, PointGe);
919 if(Intersects_AB)
920 {
921 CurrentE_Point = PointGe;
922 }
923 else
924 {
925 // No intersection found for first part of curve
926 // (CurrentA,InterMedPoint[depth]). Go to the second part
927 //
928 Second_half = true;
929 }
930 } // if did_len
931
932 if( (Second_half)&&(depth!=0) )
933 {
934 // Second part of curve (InterMed[depth],Intermed[depth-1]) )
935 // On the depth-1 level normally we are on the 'second_half'
936
937 Second_half = true;
938
939 // Find new trial intersection point needed at start of the loop
940 //
941 SubStart_PointVelocity = *ptrInterMedFT[depth];
942 CurrentA_PointVelocity = *ptrInterMedFT[depth];
943 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
944 // Ensure that the new endpoints are not further apart in space
945 // than on the curve due to different errors in the integration
946 //
947 G4double linDistSq, curveDist;
948 linDistSq = ( CurrentB_PointVelocity.GetPosition()
949 - CurrentA_PointVelocity.GetPosition() ).mag2();
950 curveDist = CurrentB_PointVelocity.GetCurveLength()
951 - CurrentA_PointVelocity.GetCurveLength();
952 if( curveDist*curveDist*(1+2*fEpsilonStep ) < linDistSq )
953 {
954 // Re-integrate to obtain a new B
955 //
956 G4FieldTrack newEndPointFT=
957 ReEstimateEndpoint( CurrentA_PointVelocity,
958 CurrentB_PointVelocity,
959 linDistSq, // to avoid recalculation
960 curveDist );
961 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
962 CurrentB_PointVelocity = newEndPointFT;
963 maxNumberOfCallsToReIntegration_depth= maxNumberOfCallsToReIntegration_depth+1;
964 }
965
966 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
967 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
968 fNavigator->LocateGlobalPointWithinVolume(Point_A);
969 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
970 stepLengthAB, PointGe);
971 if(Intersects_AB)
972 {
973 CurrentE_Point = PointGe;
974 }
975 else
976 {
977 final_section = true;
978 }
979 depth--;
980 }
981 } // if(!found_aproximate_intersection)
982
983 } while ( ( ! found_approximate_intersection )
984 && ( ! there_is_no_intersection )
985 && ( substep_no <= max_substeps) ); // UNTIL found or failed
986
987 if( substep_no > max_no_seen )
988 {
989 max_no_seen = substep_no;
990 if( max_no_seen > warn_substeps )
991 {
992 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
993 }
994 }
995
996 if( ( substep_no >= max_substeps)
997 && !there_is_no_intersection
998 && !found_approximate_intersection )
999 {
1000 G4cerr << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
1001 << G4endl
1002 << " Convergence is requiring too many substeps: "
1003 << substep_no << G4endl;
1004 G4cerr << " Abandoning effort to intersect. " << G4endl;
1005 G4cerr << " Information on start & current step follows in cout."
1006 << G4endl;
1007 G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
1008 << G4endl
1009 << " Convergence is requiring too many substeps: "
1010 << substep_no << G4endl;
1011 G4cout << " Found intersection = "
1012 << found_approximate_intersection << G4endl
1013 << " Intersection exists = "
1014 << !there_is_no_intersection << G4endl;
1015 G4cout << " Start and Endpoint of Requested Step:" << G4endl;
1016 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
1017 -1.0, NewSafety, 0, 0);
1018 G4cout << G4endl;
1019 G4cout << " 'Bracketing' starting and endpoint of current Sub-Step"
1020 << G4endl;
1021 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
1022 -1.0, NewSafety, substep_no-1, 0);
1023 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
1024 -1.0, NewSafety, substep_no, 0);
1025 G4cout << G4endl;
1026
1027#ifdef FUTURE_CORRECTION
1028 // Attempt to correct the results of the method // FIX - TODO
1029
1030 if ( ! found_approximate_intersection )
1031 {
1032 recalculatedEndPoint = true;
1033 // Return the further valid intersection point -- potentially A ??
1034 // JA/19 Jan 2006
1035 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
1036
1037 G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
1038 << G4endl
1039 << " Did not convergence after " << substep_no
1040 << " substeps." << G4endl;
1041 G4cout << " The endpoint was adjused to pointA resulting"
1042 << G4endl
1043 << " from the last substep: " << CurrentA_PointVelocity
1044 << G4endl;
1045 }
1046#endif
1047
1048 G4cout.precision( 10 );
1049 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
1050 G4double full_len = CurveEndPointVelocity.GetCurveLength();
1051 G4cout << "ERROR - G4PropagatorInField::LocateIntersectionPoint()"
1052 << G4endl
1053 << " Undertaken only length: " << done_len
1054 << " out of " << full_len << " required." << G4endl;
1055 G4cout << " Remaining length = " << full_len - done_len << G4endl;
1056
1057 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
1058 "UnableToLocateIntersection", FatalException,
1059 "Too many substeps while trying to locate intersection.");
1060 }
1061 else if( substep_no >= warn_substeps )
1062 {
1063 int oldprc= G4cout.precision( 10 );
1064 G4cout << "WARNING - G4PropagatorInField::LocateIntersectionPoint()"
1065 << G4endl
1066 << " Undertaken length: "
1067 << CurrentB_PointVelocity.GetCurveLength();
1068 G4cout << " - Needed: " << substep_no << " substeps." << G4endl
1069 << " Warning level = " << warn_substeps
1070 << " and maximum substeps = " << max_substeps << G4endl;
1071 G4Exception("G4PropagatorInField::LocateIntersectionPoint()",
1072 "DifficultyToLocateIntersection", JustWarning,
1073 "Many substeps while trying to locate intersection.");
1074 G4cout.precision( oldprc );
1075 }
1076
1077 return !there_is_no_intersection; // Success or failure
1078}
1079
1080///////////////////////////////////////////////////////////////////////////
1081//
1082// Dumps status of propagator.
1083
1084void
1085G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
1086 const G4FieldTrack& CurrentFT,
1087 G4double requestStep,
1088 G4double safety,
1089 G4int stepNo,
1090 G4VPhysicalVolume* startVolume)
1091{
1092 const G4int verboseLevel= fVerboseLevel;
1093 const G4ThreeVector StartPosition = StartFT.GetPosition();
1094 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
1095 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
1096 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
1097
1098 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
1099
1100 if( ((stepNo == 0) && (verboseLevel <3))
1101 || (verboseLevel >= 3) )
1102 {
1103 static G4int noPrecision= 4;
1104 G4cout.precision(noPrecision);
1105 // G4cout.setf(ios_base::fixed,ios_base::floatfield);
1106 G4cout << std::setw( 6) << " "
1107 << std::setw( 25) << " Current Position and Direction" << " "
1108 << G4endl;
1109 G4cout << std::setw( 5) << "Step#"
1110 << std::setw(10) << " s " << " "
1111 << std::setw(10) << "X(mm)" << " "
1112 << std::setw(10) << "Y(mm)" << " "
1113 << std::setw(10) << "Z(mm)" << " "
1114 << std::setw( 7) << " N_x " << " "
1115 << std::setw( 7) << " N_y " << " "
1116 << std::setw( 7) << " N_z " << " " ;
1117 // << G4endl;
1118 G4cout // << " >>> "
1119 << std::setw( 7) << " Delta|N|" << " "
1120 // << std::setw( 7) << " Delta(N_z) " << " "
1121 << std::setw( 9) << "StepLen" << " "
1122 << std::setw(12) << "StartSafety" << " "
1123 << std::setw( 9) << "PhsStep" << " ";
1124 if( startVolume ) {
1125 G4cout << std::setw(18) << "NextVolume" << " ";
1126 }
1127 G4cout << G4endl;
1128 }
1129 if((stepNo == 0) && (verboseLevel <=3)){
1130 // Recurse to print the start values
1131 //
1132 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
1133 }
1134 if( verboseLevel <= 3 )
1135 {
1136 if( stepNo >= 0)
1137 G4cout << std::setw( 4) << stepNo << " ";
1138 else
1139 G4cout << std::setw( 5) << "Start" ;
1140 G4cout.precision(8);
1141 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
1142 G4cout.precision(8);
1143 G4cout << std::setw(10) << CurrentPosition.x() << " "
1144 << std::setw(10) << CurrentPosition.y() << " "
1145 << std::setw(10) << CurrentPosition.z() << " ";
1146 G4cout.precision(4);
1147 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
1148 << std::setw( 7) << CurrentUnitVelocity.y() << " "
1149 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
1150 // G4cout << G4endl;
1151 // G4cout << " >>> " ;
1152 G4cout.precision(3);
1153 G4cout << std::setw( 7) << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag() << " ";
1154 // << std::setw( 7) << CurrentUnitVelocity.z() - InitialUnitVelocity.z() << " ";
1155 G4cout << std::setw( 9) << step_len << " ";
1156 G4cout << std::setw(12) << safety << " ";
1157 if( requestStep != -1.0 )
1158 G4cout << std::setw( 9) << requestStep << " ";
1159 else
1160 G4cout << std::setw( 9) << "Init/NotKnown" << " ";
1161
1162 if( startVolume != 0)
1163 {
1164 G4cout << std::setw(12) << startVolume->GetName() << " ";
1165 }
1166#if 0
1167 else
1168 {
1169 if( step_len != -1 )
1170 G4cout << std::setw(12) << "OutOfWorld" << " ";
1171 else
1172 G4cout << std::setw(12) << "NotGiven" << " ";
1173 }
1174#endif
1175
1176 G4cout << G4endl;
1177 }
1178 else // if( verboseLevel > 3 )
1179 {
1180 // Multi-line output
1181
1182 G4cout << "Step taken was " << step_len
1183 << " out of PhysicalStep= " << requestStep << G4endl;
1184 G4cout << "Final safety is: " << safety << G4endl;
1185
1186 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
1187 << G4endl;
1188 G4cout << G4endl;
1189 }
1190}
1191
1192///////////////////////////////////////////////////////////////////////////
1193//
1194// Prints Step diagnostics
1195
1196void
1197G4PropagatorInField::PrintStepLengthDiagnostic(
1198 G4double CurrentProposedStepLength,
1199 G4double decreaseFactor,
1200 G4double stepTrial,
1201 const G4FieldTrack& )
1202{
1203 G4cout << " PiF: NoZeroStep= " << fNoZeroStep
1204 << " CurrentProposedStepLength= " << CurrentProposedStepLength
1205 << " Full_curvelen_last=" << fFull_CurveLen_of_LastAttempt
1206 << " last proposed step-length= " << fLast_ProposedStepLength
1207 << " decreate factor = " << decreaseFactor
1208 << " step trial = " << stepTrial
1209 << G4endl;
1210}
1211
1212G4bool
1213G4PropagatorInField::IntersectChord( G4ThreeVector StartPointA,
1214 G4ThreeVector EndPointB,
1215 G4double &NewSafety,
1216 G4double &LinearStepLength,
1217 G4ThreeVector &IntersectionPoint
1218 )
1219{
1220 // Calculate the direction and length of the chord AB
1221 G4ThreeVector ChordAB_Vector = EndPointB - StartPointA;
1222 G4double ChordAB_Length = ChordAB_Vector.mag(); // Magnitude (norm)
1223 G4ThreeVector ChordAB_Dir = ChordAB_Vector.unit();
1224 G4bool intersects;
1225
1226 G4ThreeVector OriginShift = StartPointA - fPreviousSftOrigin ;
1227 G4double MagSqShift = OriginShift.mag2() ;
1228 G4double currentSafety;
1229 G4bool doCallNav= false;
1230
1231 if( MagSqShift >= sqr(fPreviousSafety) )
1232 {
1233 currentSafety = 0.0 ;
1234 }else{
1235 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
1236 }
1237
1238 if( fUseSafetyForOptimisation && (ChordAB_Length <= currentSafety) )
1239 {
1240 // The Step is guaranteed to be taken
1241
1242 LinearStepLength = ChordAB_Length;
1243 intersects = false;
1244
1245 NewSafety= currentSafety;
1246
1247#if 0
1248 G4cout << " G4PropagatorInField does not call Navigator::ComputeStep " << G4endl ;
1249 G4cout << " step= " << LinearStepLength << " safety= " << NewSafety << G4endl;
1250 G4cout << " safety: Origin = " << fPreviousSftOrigin << " val= " << fPreviousSafety << G4endl;
1251#endif
1252 }
1253 else
1254 {
1255 doCallNav= true;
1256 // Check whether any volumes are encountered by the chord AB
1257
1258 // G4cout << " G4PropagatorInField calling Navigator::ComputeStep " << G4endl ;
1259
1260 LinearStepLength =
1261 fNavigator->ComputeStep( StartPointA, ChordAB_Dir,
1262 ChordAB_Length, NewSafety );
1263 intersects = (LinearStepLength <= ChordAB_Length);
1264 // G4Navigator contracts to return k_infinity if len==asked
1265 // and it did not find a surface boundary at that length
1266 LinearStepLength = std::min( LinearStepLength, ChordAB_Length);
1267
1268 // G4cout << " G4PiF got step= " << LinearStepLength << " safety= " << NewSafety << G4endl;
1269
1270 // Save the last calculated safety!
1271 fPreviousSftOrigin = StartPointA;
1272 fPreviousSafety= NewSafety;
1273
1274 if( intersects ){
1275 // Intersection Point of chord AB and either volume A's surface
1276 // or a daughter volume's surface ..
1277 IntersectionPoint = StartPointA + LinearStepLength * ChordAB_Dir;
1278 }
1279 }
1280
1281#ifdef DEBUG_INTERSECTS_CHORD
1282 // printIntersection(
1283 // StartPointA, EndPointB, LinearStepLength, IntersectionPoint, NewSafety
1284
1285 G4cout << " G4PropagatorInField::IntersectChord reports " << G4endl;
1286 G4cout << " PiF-IC> "
1287 << "Start=" << std::setw(12) << StartPointA << " "
1288 << "End= " << std::setw(8) << EndPointB << " "
1289 << "StepIn=" << std::setw(8) << LinearStepLength << " "
1290 << "NewSft=" << std::setw(8) << NewSafety << " "
1291 << "CallNav=" << doCallNav << " "
1292 << "Intersects " << intersects << " ";
1293 if( intersects )
1294 G4cout << "IntrPt=" << std::setw(8) << IntersectionPoint << " " ;
1295 G4cout << G4endl;
1296#endif
1297
1298 return intersects;
1299}
1300
1301// --------------------- oooo000000000000oooo ----------------------------
1302
1303G4FieldTrack G4PropagatorInField::
1304ReEstimateEndpoint( const G4FieldTrack &CurrentStateA,
1305 const G4FieldTrack &EstimatedEndStateB,
1306 G4double linearDistSq,
1307 G4double curveDist
1308 )
1309{
1310 // G4double checkCurveDist= EstimatedEndStateB.GetCurveLength()
1311 // - CurrentStateA.GetCurveLength();
1312 // G4double checkLinDistSq= (EstimatedEndStateB.GetPosition()
1313 // - CurrentStateA.GetPosition() ).mag2();
1314
1315 G4FieldTrack newEndPoint( CurrentStateA );
1316 G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
1317
1318 G4FieldTrack retEndPoint( CurrentStateA );
1319 G4bool goodAdvance;
1320 G4int itrial=0;
1321 const G4int no_trials= 20;
1322
1323 G4double endCurveLen= EstimatedEndStateB.GetCurveLength();
1324 do
1325 {
1326 G4double currentCurveLen= newEndPoint.GetCurveLength();
1327 G4double advanceLength= endCurveLen - currentCurveLen ;
1328 if (std::abs(advanceLength)<kCarTolerance)
1329 {
1330 advanceLength=(EstimatedEndStateB.GetPosition()
1331 -newEndPoint.GetPosition()).mag();
1332 }
1333 goodAdvance=
1334 integrDriver->AccurateAdvance(newEndPoint, advanceLength, fEpsilonStep);
1335 // ***************
1336 }
1337 while( !goodAdvance && (++itrial < no_trials) );
1338
1339 if( goodAdvance )
1340 {
1341 retEndPoint= newEndPoint;
1342 }
1343 else
1344 {
1345 retEndPoint= EstimatedEndStateB; // Could not improve without major work !!
1346 }
1347
1348 // All the work is done
1349 // below are some diagnostics only -- before the return!
1350 //
1351 static const G4String MethodName("G4PropagatorInField::ReEstimateEndpoint");
1352
1353#ifdef G4VERBOSE
1354 G4int latest_good_trials=0;
1355 if( itrial > 1)
1356 {
1357 if( fVerboseLevel > 0 )
1358 {
1359 G4cout << MethodName << " called - goodAdv= " << goodAdvance
1360 << " trials = " << itrial
1361 << " previous good= " << latest_good_trials
1362 << G4endl;
1363 }
1364 latest_good_trials=0;
1365 }
1366 else
1367 {
1368 latest_good_trials++;
1369 }
1370#endif
1371
1372#ifdef G4DEBUG_FIELD
1373 G4double lengthDone = newEndPoint.GetCurveLength()
1374 - CurrentStateA.GetCurveLength();
1375 if( !goodAdvance )
1376 {
1377 if( fVerboseLevel >= 3 )
1378 {
1379 G4cout << MethodName << "> AccurateAdvance failed " ;
1380 G4cout << " in " << itrial << " integration trials/steps. " << G4endl;
1381 G4cout << " It went only " << lengthDone << " instead of " << curveDist
1382 << " -- a difference of " << curveDist - lengthDone << G4endl;
1383 G4cout << " ReEstimateEndpoint> Reset endPoint to original value!"
1384 << G4endl;
1385 }
1386 }
1387
1388 static G4int noInaccuracyWarnings = 0;
1389 G4int maxNoWarnings = 10;
1390 if ( (noInaccuracyWarnings < maxNoWarnings )
1391 || (fVerboseLevel > 1) )
1392 {
1393 G4cerr << "G4PropagatorInField::LocateIntersectionPoint():"
1394 << G4endl
1395 << " Warning: Integration inaccuracy requires"
1396 << " an adjustment in the step's endpoint." << G4endl
1397 << " Two mid-points are further apart than their"
1398 << " curve length difference" << G4endl
1399 << " Dist = " << std::sqrt(linearDistSq)
1400 << " curve length = " << curveDist << G4endl;
1401 G4cerr << " Correction applied is "
1402 << (newEndPoint.GetPosition()-EstimatedEndStateB.GetPosition()).mag()
1403 << G4endl;
1404 }
1405#else
1406 // Statistics on the RMS value of the corrections
1407
1408 static G4int noCorrections=0;
1409 static G4double sumCorrectionsSq = 0;
1410 noCorrections++;
1411 if( goodAdvance )
1412 {
1413 sumCorrectionsSq += (EstimatedEndStateB.GetPosition() -
1414 newEndPoint.GetPosition()).mag2();
1415 }
1416 linearDistSq -= curveDist; // To use linearDistSq ... !
1417#endif
1418
1419 return retEndPoint;
1420}
1421
1422// Access the points which have passed through the filter. The
1423// points are stored as ThreeVectors for the initial impelmentation
1424// only (jacek 30/10/2002)
1425// Responsibility for deleting the points lies with
1426// SmoothTrajectoryPoint, which is the points' final
1427// destination. The points pointer is set to NULL, to ensure that
1428// the points are not re-used in subsequent steps, therefore THIS
1429// METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
1430
1431std::vector<G4ThreeVector>*
1432G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
1433{
1434 // NB, GimmeThePointsAndForgetThem really forgets them, so it can
1435 // only be called (exactly) once for each step.
1436
1437 if (fpTrajectoryFilter)
1438 {
1439 return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
1440 }
1441 else
1442 {
1443 return 0;
1444 }
1445}
1446
1447void
1448G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
1449{
1450 fpTrajectoryFilter = filter;
1451}
1452
1453void G4PropagatorInField::ClearPropagatorState()
1454{
1455 // Goal: Clear all memory of previous steps, cached information
1456
1457 fParticleIsLooping= false;
1458 fNoZeroStep= 0;
1459
1460 End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
1461 G4ThreeVector(0.,0.,0.),
1462 0.0,0.0,0.0,0.0,0.0);
1463 fFull_CurveLen_of_LastAttempt = -1;
1464 fLast_ProposedStepLength = -1;
1465
1466 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
1467 fPreviousSafety= 0.0;
1468}
1469
1470G4FieldManager* G4PropagatorInField::
1471FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
1472{
1473 G4FieldManager* currentFieldMgr;
1474
1475 currentFieldMgr = fDetectorFieldMgr;
1476 if( pCurrentPhysicalVolume)
1477 {
1478 G4FieldManager *newFieldMgr = 0;
1479 newFieldMgr= pCurrentPhysicalVolume->GetLogicalVolume()->GetFieldManager();
1480 if ( newFieldMgr )
1481 currentFieldMgr = newFieldMgr;
1482 }
1483 fCurrentFieldMgr= currentFieldMgr;
1484
1485 // Flag that field manager has been set.
1486 fSetFieldMgr= true;
1487
1488 return currentFieldMgr;
1489}
1490
1491G4int G4PropagatorInField::SetVerboseLevel( G4int level )
1492{
1493 G4int oldval= fVerboseLevel;
1494 fVerboseLevel= level;
1495
1496 // Forward the verbose level 'reduced' to ChordFinder,
1497 // MagIntegratorDriver ... ?
1498 //
1499 G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver();
1500 integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
1501 G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
1502
1503 return oldval;
1504}
Note: See TracBrowser for help on using the repository browser.