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

Last change on this file since 833 was 831, checked in by garnier, 17 years ago

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