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