| Geant4 Cross Reference |
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 // class G4PropagatorInField Implementation
27 //
28 // This class implements an algorithm to track a particle in a
29 // non-uniform magnetic field. It utilises an ODE solver (with
30 // the Runge - Kutta method) to evolve the particle, and drives it
31 // until the particle has traveled a set distance or it enters a new
32 // volume.
33 //
34 // 14.10.96 John Apostolakis, design and implementation
35 // 17.03.97 John Apostolakis, renaming new set functions being added
36 // ---------------------------------------------------------------------------
37
38 #include <iomanip>
39
40 #include "G4PropagatorInField.hh"
41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4ThreeVector.hh"
44 #include "G4Material.hh"
45 #include "G4VPhysicalVolume.hh"
46 #include "G4Navigator.hh"
47 #include "G4GeometryTolerance.hh"
48 #include "G4VCurvedTrajectoryFilter.hh"
49 #include "G4ChordFinder.hh"
50 #include "G4MultiLevelLocator.hh"
51
52
53 // ---------------------------------------------------------------------------
54 // Constructors and destructor
55 //
56 G4PropagatorInField::G4PropagatorInField( G4Navigator* theNavigator,
57 G4FieldManager* detectorFieldMgr,
58 G4VIntersectionLocator* vLocator )
59 : fDetectorFieldMgr(detectorFieldMgr),
60 fNavigator(theNavigator),
61 fCurrentFieldMgr(detectorFieldMgr),
62 End_PointAndTangent(G4ThreeVector(0.,0.,0.),
63 G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0)
64 {
65 fEpsilonStep = (fDetectorFieldMgr != nullptr)
66 ? fDetectorFieldMgr->GetMaximumEpsilonStep() : 1.0e-5;
67
68
69 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.);
70 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
71 fZeroStepThreshold = std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer );
72
73 fLargestAcceptableStep = 100.0 * meter; // Reduced from 1000.0 * meter
74 fMaxStepSizeMultiplier= 0.1 ; // 0.1 in git (larger for tests.) // Reduced from 100;
75 fMinBigDistance= 100. * CLHEP::mm;
76 #ifdef G4DEBUG_FIELD
77 G4cout << " PiF: Zero Step Threshold set to "
78 << fZeroStepThreshold / millimeter
79 << " mm." << G4endl;
80 G4cout << " PiF: Value of kCarTolerance = "
81 << kCarTolerance / millimeter
82 << " mm. " << G4endl;
83 fVerboseLevel = 2;
84 fVerbTracePiF = true;
85 #endif
86
87 // Defining Intersection Locator and his parameters
88 if ( vLocator == nullptr )
89 {
90 fIntersectionLocator = new G4MultiLevelLocator(theNavigator);
91 fAllocatedLocator = true;
92 }
93 else
94 {
95 fIntersectionLocator = vLocator;
96 fAllocatedLocator = false;
97 }
98 RefreshIntersectionLocator(); // Copy all relevant parameters
99 }
100
101 // ---------------------------------------------------------------------------
102 //
103 G4PropagatorInField::~G4PropagatorInField()
104 {
105 if(fAllocatedLocator) { delete fIntersectionLocator; }
106 }
107
108 // ---------------------------------------------------------------------------
109 // Update the IntersectionLocator with current parameters
110 //
111 void G4PropagatorInField::RefreshIntersectionLocator()
112 {
113 fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
114 fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
115 fIntersectionLocator->SetChordFinderFor(GetChordFinder());
116 fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
117 }
118
119 // ---------------------------------------------------------------------------
120 // Compute the next geometric Step
121 //
122 G4double G4PropagatorInField::ComputeStep(
123 G4FieldTrack& pFieldTrack,
124 G4double CurrentProposedStepLength,
125 G4double& currentSafety, // IN/OUT
126 G4VPhysicalVolume* pPhysVol,
127 G4bool canRelaxDeltaChord)
128 {
129 GetChordFinder()->OnComputeStep(&pFieldTrack);
130 const G4double deltaChord = GetChordFinder()->GetDeltaChord();
131
132 // If CurrentProposedStepLength is too small for finding Chords
133 // then return with no action (for now - TODO: some action)
134 //
135 const char* methodName = "G4PropagatorInField::ComputeStep";
136 if (CurrentProposedStepLength<kCarTolerance)
137 {
138 return kInfinity;
139 }
140
141 // Introducing smooth trajectory display (jacek 01/11/2002)
142 //
143 if (fpTrajectoryFilter != nullptr)
144 {
145 fpTrajectoryFilter->CreateNewTrajectorySegment();
146 }
147
148 fFirstStepInVolume = fNewTrack ? true : fLastStepInVolume;
149 fLastStepInVolume = false;
150 fNewTrack = false;
151
152 if( fVerboseLevel > 2 )
153 {
154 G4cout << methodName << " called" << G4endl;
155 G4cout << " Starting FT: " << pFieldTrack;
156 G4cout << " Requested length = " << CurrentProposedStepLength << G4endl;
157 G4cout << " PhysVol = ";
158 if( pPhysVol != nullptr )
159 {
160 G4cout << pPhysVol->GetName() << G4endl;
161 }
162 else
163 {
164 G4cout << " N/A ";
165 }
166 G4cout << G4endl;
167 }
168
169 // Parameters for adaptive Runge-Kutta integration
170
171 G4double h_TrialStepSize; // 1st Step Size
172 G4double TruePathLength = CurrentProposedStepLength;
173 G4double StepTaken = 0.0;
174 G4double s_length_taken, epsilon;
175 G4bool intersects;
176 G4bool first_substep = true;
177
178 G4double NewSafety;
179 fParticleIsLooping = false;
180
181 // If not yet done,
182 // Set the field manager to the local one if the volume has one,
183 // or to the global one if not
184 //
185 if( !fSetFieldMgr )
186 {
187 fCurrentFieldMgr = FindAndSetFieldManager( pPhysVol );
188 }
189 fSetFieldMgr = false; // For next call, the field manager must be set again
190
191 G4FieldTrack CurrentState(pFieldTrack);
192 G4FieldTrack OriginalState = CurrentState;
193
194 // If the Step length is "infinite", then an approximate-maximum Step
195 // length (used to calculate the relative accuracy) must be guessed
196 //
197 if( CurrentProposedStepLength >= fLargestAcceptableStep )
198 {
199 G4ThreeVector StartPointA, VelocityUnit;
200 StartPointA = pFieldTrack.GetPosition();
201 VelocityUnit = pFieldTrack.GetMomentumDir();
202
203 G4double trialProposedStep = fMaxStepSizeMultiplier * ( fMinBigDistance +
204 fNavigator->GetWorldVolume()->GetLogicalVolume()->
205 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
206 CurrentProposedStepLength = std::min( trialProposedStep,
207 fLargestAcceptableStep );
208 }
209 epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
210 G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
211 G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();
212 if( epsilon < epsilonMin ) { epsilon = epsilonMin; }
213 if( epsilon > epsilonMax ) { epsilon = epsilonMax; }
214 SetEpsilonStep( epsilon );
215
216 // Values for Intersection Locator has to be updated on each call for the
217 // case that CurrentFieldManager has changed from the one of previous step
218 //
219 RefreshIntersectionLocator();
220
221 // Shorten the proposed step in case of earlier problems (zero steps)
222 //
223 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
224 {
225 G4double stepTrial;
226
227 stepTrial = fFull_CurveLen_of_LastAttempt;
228 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
229 {
230 stepTrial = fLast_ProposedStepLength;
231 }
232
233 G4double decreaseFactor = 0.9; // Unused default
234 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
235 && (stepTrial > 100.0*fZeroStepThreshold) )
236 {
237 // Attempt quick convergence
238 //
239 decreaseFactor= 0.25;
240 }
241 else
242 {
243 // We are in significant difficulties, probably at a boundary that
244 // is either geometrically sharp or between very different materials.
245 // Careful decreases to cope with tolerance are required
246 //
247 if( stepTrial > 100.0*fZeroStepThreshold ) {
248 decreaseFactor = 0.35; // Try decreasing slower
249 } else if( stepTrial > 30.0*fZeroStepThreshold ) {
250 decreaseFactor= 0.5; // Try yet slower decrease
251 } else if( stepTrial > 10.0*fZeroStepThreshold ) {
252 decreaseFactor= 0.75; // Try even slower decreases
253 } else {
254 decreaseFactor= 0.9; // Try very slow decreases
255 }
256 }
257 stepTrial *= decreaseFactor;
258
259 #ifdef G4DEBUG_FIELD
260 if( fVerboseLevel > 2
261 || (fNoZeroStep >= fSevereActionThreshold_NoZeroSteps) )
262 {
263 G4cerr << " " << methodName
264 << " Decreasing step after " << fNoZeroStep << " zero steps "
265 << " - in volume " << pPhysVol;
266 if( pPhysVol )
267 G4cerr << " with name " << pPhysVol->GetName();
268 else
269 G4cerr << " i.e. *unknown* volume.";
270 G4cerr << G4endl;
271 PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
272 stepTrial, pFieldTrack);
273 }
274 #endif
275 if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
276 {
277 std::ostringstream message;
278 message << "Particle abandoned due to lack of progress in field."
279 << G4endl
280 << " Properties : " << pFieldTrack << G4endl
281 << " Attempting a zero step = " << stepTrial << G4endl
282 << " while attempting to progress after " << fNoZeroStep
283 << " trial steps. Will abandon step.";
284 G4Exception(methodName, "GeomNav1002", JustWarning, message);
285 fParticleIsLooping = true;
286 return 0; // = stepTrial;
287 }
288 if( stepTrial < CurrentProposedStepLength )
289 {
290 CurrentProposedStepLength = stepTrial;
291 }
292 }
293 fLast_ProposedStepLength = CurrentProposedStepLength;
294
295 G4int do_loop_count = 0;
296 do // Loop checking, 07.10.2016, JA
297 {
298 G4FieldTrack SubStepStartState = CurrentState;
299 G4ThreeVector SubStartPoint = CurrentState.GetPosition();
300
301 if(!first_substep)
302 {
303 if( fVerboseLevel > 4 )
304 {
305 G4cout << " PiF: Calling Nav/Locate Global Point within-Volume "
306 << G4endl;
307 }
308 fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
309 }
310
311 // How far to attempt to move the particle !
312 //
313 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
314
315 if (canRelaxDeltaChord &&
316 fIncreaseChordDistanceThreshold > 0 &&
317 do_loop_count > fIncreaseChordDistanceThreshold &&
318 do_loop_count % fIncreaseChordDistanceThreshold == 0)
319 {
320 GetChordFinder()->SetDeltaChord(
321 GetChordFinder()->GetDeltaChord() * 2.0
322 );
323 }
324
325 // Integrate as far as "chord miss" rule allows.
326 //
327 s_length_taken = GetChordFinder()->AdvanceChordLimited(
328 CurrentState, // Position & velocity
329 h_TrialStepSize,
330 fEpsilonStep,
331 fPreviousSftOrigin,
332 fPreviousSafety );
333 // CurrentState is now updated with the final position and velocity
334
335 fFull_CurveLen_of_LastAttempt = s_length_taken;
336
337 G4ThreeVector EndPointB = CurrentState.GetPosition();
338 G4ThreeVector InterSectionPointE;
339 G4double LinearStepLength;
340
341 // Intersect chord AB with geometry
342 //
343 intersects= IntersectChord( SubStartPoint, EndPointB,
344 NewSafety, LinearStepLength,
345 InterSectionPointE );
346 // E <- Intersection Point of chord AB and either volume A's surface
347 // or a daughter volume's surface ..
348
349 if( first_substep )
350 {
351 currentSafety = NewSafety;
352 } // Updating safety in other steps is potential future extention
353
354 if( intersects )
355 {
356 G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
357
358 // Find the intersection point of AB true path with the surface
359 // of vol(A), if it exists. Start with point E as first "estimate".
360 G4bool recalculatedEndPt = false;
361
362 G4bool found_intersection = fIntersectionLocator->
363 EstimateIntersectionPoint( SubStepStartState, CurrentState,
364 InterSectionPointE, IntersectPointVelct_G,
365 recalculatedEndPt, fPreviousSafety,
366 fPreviousSftOrigin);
367 intersects = found_intersection;
368 if( found_intersection )
369 {
370 End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
371 StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
372 - OriginalState.GetCurveLength();
373 }
374 else
375 {
376 // Either "minor" chords do not intersect
377 // or else stopped (due to too many steps)
378 //
379 if( recalculatedEndPt )
380 {
381 G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
382 G4double endExpected = CurrentState.GetCurveLength();
383
384 // Detect failure - due to too many steps
385 G4bool shortEnd = endAchieved
386 < (endExpected*(1.0-CLHEP::perMillion));
387
388 G4double stepAchieved = endAchieved
389 - SubStepStartState.GetCurveLength();
390
391 // Update remaining state - must work for 'full' step or
392 // abandonned intersection
393 //
394 CurrentState = IntersectPointVelct_G;
395 s_length_taken = stepAchieved;
396 if( shortEnd )
397 {
398 fParticleIsLooping = true;
399 }
400 }
401 }
402 }
403 if( !intersects )
404 {
405 StepTaken += s_length_taken;
406
407 if (fpTrajectoryFilter != nullptr) // For smooth trajectory display (jacek 1/11/2002)
408 {
409 fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
410 }
411 }
412 first_substep = false;
413
414 #ifdef G4DEBUG_FIELD
415 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
416 {
417 if( fNoZeroStep > fSevereActionThreshold_NoZeroSteps )
418 G4cout << " Above 'Severe Action' threshold -- for Zero steps. ";
419 else
420 G4cout << " Above 'action' threshold -- for Zero steps. ";
421 G4cout << " Number of zero steps = " << fNoZeroStep << G4endl;
422 printStatus( SubStepStartState, // or OriginalState,
423 CurrentState, CurrentProposedStepLength,
424 NewSafety, do_loop_count, pPhysVol );
425 }
426 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
427 {
428 if( do_loop_count == fMax_loop_count-9 )
429 {
430 G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
431 << " Difficult track - taking many sub steps." << G4endl;
432 printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength,
433 NewSafety, 0, pPhysVol );
434 }
435 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
436 NewSafety, do_loop_count, pPhysVol );
437 }
438 #endif
439
440 ++do_loop_count;
441
442 } while( (!intersects )
443 && (!fParticleIsLooping)
444 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
445 && ( do_loop_count < fMax_loop_count ) );
446
447 if( do_loop_count >= fMax_loop_count
448 && (StepTaken + kCarTolerance < CurrentProposedStepLength) )
449 {
450 fParticleIsLooping = true;
451 }
452 if ( ( fParticleIsLooping ) && (fVerboseLevel > 0) )
453 {
454 ReportLoopingParticle( do_loop_count, StepTaken,
455 CurrentProposedStepLength, methodName,
456 CurrentState.GetMomentum(), pPhysVol );
457 }
458
459 if( !intersects )
460 {
461 // Chord AB or "minor chords" do not intersect
462 // B is the endpoint Step of the current Step.
463 //
464 End_PointAndTangent = CurrentState;
465 TruePathLength = StepTaken; // Original code
466
467 // Tried the following to avoid potential issue with round-off error
468 // - but has issues... Suppressing this change JA 2015/05/02
469 // TruePathLength = CurrentProposedStepLength;
470 }
471 fLastStepInVolume = intersects;
472
473 // Set pFieldTrack to the return value
474 //
475 pFieldTrack = End_PointAndTangent;
476
477 #ifdef G4VERBOSE
478 // Check that "s" is correct
479 //
480 if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
481 - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
482 {
483 std::ostringstream message;
484 message << "Curve length mis-match between original state "
485 << "and proposed endpoint of propagation." << G4endl
486 << " The curve length of the endpoint should be: "
487 << OriginalState.GetCurveLength() + TruePathLength << G4endl
488 << " and it is instead: "
489 << End_PointAndTangent.GetCurveLength() << "." << G4endl
490 << " A difference of: "
491 << OriginalState.GetCurveLength() + TruePathLength
492 - End_PointAndTangent.GetCurveLength() << G4endl
493 << " Original state = " << OriginalState << G4endl
494 << " Proposed state = " << End_PointAndTangent;
495 G4Exception(methodName, "GeomNav0003", FatalException, message);
496 }
497 #endif
498
499 if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
500 {
501 fNoZeroStep = 0;
502 }
503 else
504 {
505 // In particular anomalous cases, we can get repeated zero steps
506 // We identify these cases and take corrective action when they occur.
507 //
508 if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
509 {
510 ++fNoZeroStep;
511 }
512 else
513 {
514 fNoZeroStep = 0;
515 }
516 }
517 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
518 {
519 fParticleIsLooping = true;
520 ReportStuckParticle( fNoZeroStep, CurrentProposedStepLength,
521 fFull_CurveLen_of_LastAttempt, pPhysVol );
522 fNoZeroStep = 0;
523 }
524
525 GetChordFinder()->SetDeltaChord(deltaChord);
526 return TruePathLength;
527 }
528
529 // ---------------------------------------------------------------------------
530 // Dumps status of propagator
531 //
532 void
533 G4PropagatorInField::printStatus( const G4FieldTrack& StartFT,
534 const G4FieldTrack& CurrentFT,
535 G4double requestStep,
536 G4double safety,
537 G4int stepNo,
538 G4VPhysicalVolume* startVolume)
539 {
540 const G4int verboseLevel = fVerboseLevel;
541 const G4ThreeVector StartPosition = StartFT.GetPosition();
542 const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
543 const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
544 const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
545
546 G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
547
548 G4long oldprec; // cout/cerr precision settings
549
550 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
551 {
552 oldprec = G4cout.precision(4);
553 G4cout << std::setw( 5) << "Step#"
554 << std::setw(10) << " s " << " "
555 << std::setw(10) << "X(mm)" << " "
556 << std::setw(10) << "Y(mm)" << " "
557 << std::setw(10) << "Z(mm)" << " "
558 << std::setw( 7) << " N_x " << " "
559 << std::setw( 7) << " N_y " << " "
560 << std::setw( 7) << " N_z " << " " ;
561 G4cout << std::setw( 7) << " Delta|N|" << " "
562 << std::setw( 9) << "StepLen" << " "
563 << std::setw(12) << "StartSafety" << " "
564 << std::setw( 9) << "PhsStep" << " ";
565 if( startVolume != nullptr )
566 { G4cout << std::setw(18) << "NextVolume" << " "; }
567 G4cout.precision(oldprec);
568 G4cout << G4endl;
569 }
570 if((stepNo == 0) && (verboseLevel <=3))
571 {
572 // Recurse to print the start values
573 //
574 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
575 }
576 if( verboseLevel <= 3 )
577 {
578 if( stepNo >= 0)
579 { G4cout << std::setw( 4) << stepNo << " "; }
580 else
581 { G4cout << std::setw( 5) << "Start" ; }
582 oldprec = G4cout.precision(8);
583 G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
584 G4cout.precision(8);
585 G4cout << std::setw(10) << CurrentPosition.x() << " "
586 << std::setw(10) << CurrentPosition.y() << " "
587 << std::setw(10) << CurrentPosition.z() << " ";
588 G4cout.precision(4);
589 G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
590 << std::setw( 7) << CurrentUnitVelocity.y() << " "
591 << std::setw( 7) << CurrentUnitVelocity.z() << " ";
592 G4cout.precision(3);
593 G4cout << std::setw( 7)
594 << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
595 G4cout << std::setw( 9) << step_len << " ";
596 G4cout << std::setw(12) << safety << " ";
597 if( requestStep != -1.0 )
598 { G4cout << std::setw( 9) << requestStep << " "; }
599 else
600 { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
601 if( startVolume != nullptr)
602 { G4cout << std::setw(12) << startVolume->GetName() << " "; }
603 G4cout.precision(oldprec);
604 G4cout << G4endl;
605 }
606 else // if( verboseLevel > 3 )
607 {
608 // Multi-line output
609
610 G4cout << "Step taken was " << step_len
611 << " out of PhysicalStep = " << requestStep << G4endl;
612 G4cout << "Final safety is: " << safety << G4endl;
613 G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
614 << G4endl;
615 G4cout << G4endl;
616 }
617 }
618
619 // ---------------------------------------------------------------------------
620 // Prints Step diagnostics
621 //
622 void
623 G4PropagatorInField::PrintStepLengthDiagnostic(
624 G4double CurrentProposedStepLength,
625 G4double decreaseFactor,
626 G4double stepTrial,
627 const G4FieldTrack& )
628 {
629 G4long iprec= G4cout.precision(8);
630 G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
631 << " " << std::setw(20) << " CurrentProposed len "
632 << " " << std::setw(18) << " Full_curvelen_last"
633 << " " << std::setw(18) << " last proposed len "
634 << " " << std::setw(18) << " decrease factor "
635 << " " << std::setw(15) << " step trial "
636 << G4endl;
637
638 G4cout << " " << std::setw(10) << fNoZeroStep << " "
639 << " " << std::setw(20) << CurrentProposedStepLength
640 << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
641 << " " << std::setw(18) << fLast_ProposedStepLength
642 << " " << std::setw(18) << decreaseFactor
643 << " " << std::setw(15) << stepTrial
644 << G4endl;
645 G4cout.precision( iprec );
646 }
647
648 // Access the points which have passed through the filter. The
649 // points are stored as ThreeVectors for the initial impelmentation
650 // only (jacek 30/10/2002)
651 // Responsibility for deleting the points lies with
652 // SmoothTrajectoryPoint, which is the points' final
653 // destination. The points pointer is set to NULL, to ensure that
654 // the points are not re-used in subsequent steps, therefore THIS
655 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
656
657 std::vector<G4ThreeVector>*
658 G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
659 {
660 // NB, GimmeThePointsAndForgetThem really forgets them, so it can
661 // only be called (exactly) once for each step.
662
663 if (fpTrajectoryFilter != nullptr)
664 {
665 return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
666 }
667 return nullptr;
668 }
669
670 // ---------------------------------------------------------------------------
671 //
672 void
673 G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
674 {
675 fpTrajectoryFilter = filter;
676 }
677
678 // ---------------------------------------------------------------------------
679 //
680 void G4PropagatorInField::ClearPropagatorState()
681 {
682 // Goal: Clear all memory of previous steps, cached information
683
684 fParticleIsLooping = false;
685 fNoZeroStep = 0;
686
687 fSetFieldMgr = false; // Has field-manager been set for the current step?
688 fEpsilonStep= 1.0e-5; // Relative accuracy of current Step
689
690 End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
691 G4ThreeVector(0.,0.,0.),
692 0.0,0.0,0.0,0.0,0.0);
693 fFull_CurveLen_of_LastAttempt = -1;
694 fLast_ProposedStepLength = -1;
695
696 fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
697 fPreviousSafety= 0.0;
698
699 fNewTrack = true;
700 }
701
702 // ---------------------------------------------------------------------------
703 //
704 G4FieldManager* G4PropagatorInField::
705 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume )
706 {
707 G4FieldManager* currentFieldMgr;
708
709 currentFieldMgr = fDetectorFieldMgr;
710 if( pCurrentPhysicalVolume != nullptr )
711 {
712 G4FieldManager *pRegionFieldMgr = nullptr, *localFieldMgr = nullptr;
713 G4LogicalVolume* pLogicalVol = pCurrentPhysicalVolume->GetLogicalVolume();
714
715 if( pLogicalVol != nullptr )
716 {
717 // Value for Region, if any, overrides
718 //
719 G4Region* pRegion = pLogicalVol->GetRegion();
720 if( pRegion != nullptr )
721 {
722 pRegionFieldMgr = pRegion->GetFieldManager();
723 if( pRegionFieldMgr != nullptr )
724 {
725 currentFieldMgr= pRegionFieldMgr;
726 }
727 }
728
729 // 'Local' Value from logical volume, if any, overrides
730 //
731 localFieldMgr = pLogicalVol->GetFieldManager();
732 if ( localFieldMgr != nullptr )
733 {
734 currentFieldMgr = localFieldMgr;
735 }
736 }
737 }
738 fCurrentFieldMgr = currentFieldMgr;
739
740 // Flag that field manager has been set
741 //
742 fSetFieldMgr = true;
743
744 return currentFieldMgr;
745 }
746
747 // ---------------------------------------------------------------------------
748 //
749 G4int G4PropagatorInField::SetVerboseLevel( G4int level )
750 {
751 G4int oldval = fVerboseLevel;
752 fVerboseLevel = level;
753
754 // Forward the verbose level 'reduced' to ChordFinder,
755 // MagIntegratorDriver ... ?
756 //
757 auto integrDriver = GetChordFinder()->GetIntegrationDriver();
758 integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
759 G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
760
761 return oldval;
762 }
763
764 // ---------------------------------------------------------------------------
765 //
766 void G4PropagatorInField::ReportLoopingParticle( G4int count,
767 G4double StepTaken,
768 G4double StepRequested,
769 const char* methodName,
770 const G4ThreeVector& momentumVec,
771 G4VPhysicalVolume* pPhysVol )
772 {
773 std::ostringstream message;
774 G4double fraction = StepTaken / StepRequested;
775 message << " Unfinished integration of track (likely looping particle) "
776 << " of momentum " << momentumVec << " ( magnitude = "
777 << momentumVec.mag() << " ) " << G4endl
778 << " after " << count << " field substeps "
779 << " totaling " << std::setprecision(12) << StepTaken / mm << " mm "
780 << " out of requested step " << std::setprecision(12)
781 << StepRequested / mm << " mm ";
782 message << " a fraction of ";
783 G4int prec = 4;
784 if( fraction > 0.99 )
785 {
786 prec = 7;
787 }
788 else
789 {
790 if (fraction > 0.97 ) { prec = 5; }
791 }
792 message << std::setprecision(prec)
793 << 100. * StepTaken / StepRequested << " % " << G4endl ;
794 if( pPhysVol != nullptr )
795 {
796 message << " in volume " << pPhysVol->GetName() ;
797 auto material = pPhysVol->GetLogicalVolume()->GetMaterial();
798 if( material != nullptr )
799 {
800 message << " with material " << material->GetName()
801 << " ( density = "
802 << material->GetDensity() / ( g/(cm*cm*cm) ) << " g / cm^3 ) ";
803 }
804 }
805 else
806 {
807 message << " in unknown (null) volume. " ;
808 }
809 G4Exception(methodName, "GeomNav1002", JustWarning, message);
810 }
811
812 // ---------------------------------------------------------------------------
813 //
814 void G4PropagatorInField::ReportStuckParticle( G4int noZeroSteps,
815 G4double proposedStep,
816 G4double lastTriedStep,
817 G4VPhysicalVolume* physVol )
818 {
819 std::ostringstream message;
820 message << "Particle is stuck; it will be killed." << G4endl
821 << " Zero progress for " << noZeroSteps << " attempted steps."
822 << G4endl
823 << " Proposed Step is " << proposedStep
824 << " but Step Taken is "<< lastTriedStep << G4endl;
825 if( physVol != nullptr )
826 {
827 message << " in volume " << physVol->GetName() ;
828 }
829 else
830 {
831 message << " in unknown or null volume. " ;
832 }
833 G4Exception("G4PropagatorInField::ComputeStep()",
834 "GeomNav1002", JustWarning, message);
835 }
836
837 // ------------------------------------------------------------------------
838
839 // ----------------------------------------------
840 // Methods to alter Parameters
841 // ----------------------------------------------
842
843 // Was a data member (of an object) -- now moved to class member
844 G4double G4PropagatorInField::GetLargestAcceptableStep()
845 {
846 return fLargestAcceptableStep;
847 }
848
849 // ------------------------------------------------------------------------
850 //
851 void G4PropagatorInField::SetLargestAcceptableStep( G4double newBigDist )
852 {
853 if( fLargestAcceptableStep>0.0 )
854 {
855 fLargestAcceptableStep = newBigDist;
856 }
857 }
858
859 // ---------------------------------------------------------------------------
860
861 G4double G4PropagatorInField::GetMaxStepSizeMultiplier()
862 {
863 return fMaxStepSizeMultiplier;
864 }
865
866 // ---------------------------------------------------------------------------
867
868 void G4PropagatorInField::SetMaxStepSizeMultiplier(G4double vm)
869 {
870 fMaxStepSizeMultiplier=vm;
871 }
872
873 // ---------------------------------------------------------------------------
874
875 G4double G4PropagatorInField::GetMinBigDistance()
876 {
877 return fMinBigDistance;
878 }
879
880 // ---------------------------------------------------------------------------
881
882 void G4PropagatorInField::SetMinBigDistance(G4double val)
883 {
884 fMinBigDistance= val;
885 }
886