// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // // // This class implements an algorithm to track a particle in a // non-uniform magnetic field. It utilises an ODE solver (with // the Runge - Kutta method) to evolve the particle, and drives it // until the particle has traveled a set distance or it enters a new // volume. // // 14.10.96 John Apostolakis, design and implementation // 17.03.97 John Apostolakis, renaming new set functions being added // // --------------------------------------------------------------------------- #include "G4PropagatorInField.hh" #include "G4ios.hh" #include #include "G4ThreeVector.hh" #include "G4VPhysicalVolume.hh" #include "G4Navigator.hh" #include "G4GeometryTolerance.hh" #include "G4VCurvedTrajectoryFilter.hh" #include "G4ChordFinder.hh" #include "G4BrentLocator.hh" /////////////////////////////////////////////////////////////////////////// // // Constructors and destructor G4PropagatorInField::G4PropagatorInField( G4Navigator *theNavigator, G4FieldManager *detectorFieldMgr, G4VIntersectionLocator *vLocator ) : fDetectorFieldMgr(detectorFieldMgr), fCurrentFieldMgr(detectorFieldMgr), fNavigator(theNavigator), End_PointAndTangent(G4ThreeVector(0.,0.,0.), G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0), fParticleIsLooping(false), fVerboseLevel(0), fMax_loop_count(1000), fNoZeroStep(0), fCharge(0.0), fInitialMomentumModulus(0.0), fMass(0.0), fUseSafetyForOptimisation(true), // (false) is less sensitive to incorrect safety fSetFieldMgr(false), fpTrajectoryFilter( 0 ) { if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();} else { fEpsilonStep= 1.0e-5; } fActionThreshold_NoZeroSteps = 2; fSevereActionThreshold_NoZeroSteps = 10; fAbandonThreshold_NoZeroSteps = 50; fFull_CurveLen_of_LastAttempt = -1; fLast_ProposedStepLength = -1; fLargestAcceptableStep = 1000.0 * meter; fPreviousSftOrigin= G4ThreeVector(0.,0.,0.); fPreviousSafety= 0.0; kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); // Definding Intersection Locator and his parameters if(vLocator==0){ fIntersectionLocator= new G4BrentLocator(theNavigator); fAllocatedLocator=true; }else{ fIntersectionLocator=vLocator; fAllocatedLocator=false; } fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep); fIntersectionLocator->SetDeltaIntersectionFor(GetDeltaIntersection()); fIntersectionLocator->SetChordFinderFor(GetChordFinder()); fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation); } G4PropagatorInField::~G4PropagatorInField() { if(fAllocatedLocator)delete fIntersectionLocator; } /////////////////////////////////////////////////////////////////////////// // // Compute the next geometric Step G4double G4PropagatorInField::ComputeStep( G4FieldTrack& pFieldTrack, G4double CurrentProposedStepLength, G4double& currentSafety, // IN/OUT G4VPhysicalVolume* pPhysVol) { // If CurrentProposedStepLength is too small for finding Chords // then return with no action (for now - TODO: some action) // if(CurrentProposedStepLengthCreateNewTrajectorySegment(); } // Parameters for adaptive Runge-Kutta integration G4double h_TrialStepSize; // 1st Step Size G4double TruePathLength = CurrentProposedStepLength; G4double StepTaken = 0.0; G4double s_length_taken, epsilon ; G4bool intersects; G4bool first_substep = true; G4double NewSafety; fParticleIsLooping = false; // If not yet done, // Set the field manager to the local one if the volume has one, // or to the global one if not // if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol ); // For the next call, the field manager must again be set fSetFieldMgr= false; GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass); // Values for Intersection Locator has to be updated on each call // because the CurrentFieldManager changes fIntersectionLocator->SetChordFinderFor(GetChordFinder()); fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation); fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep); fIntersectionLocator->SetDeltaIntersectionFor(GetDeltaIntersection()); G4FieldTrack CurrentState(pFieldTrack); G4FieldTrack OriginalState = CurrentState; // If the Step length is "infinite", then an approximate-maximum Step // length (used to calculate the relative accuracy) must be guessed. // if( CurrentProposedStepLength >= fLargestAcceptableStep ) { G4ThreeVector StartPointA, VelocityUnit; StartPointA = pFieldTrack.GetPosition(); VelocityUnit = pFieldTrack.GetMomentumDir(); G4double trialProposedStep = 1.e2 * ( 10.0 * cm + fNavigator->GetWorldVolume()->GetLogicalVolume()-> GetSolid()->DistanceToOut(StartPointA, VelocityUnit) ); CurrentProposedStepLength= std::min( trialProposedStep, fLargestAcceptableStep ); } epsilon = GetDeltaOneStep() / CurrentProposedStepLength; // G4double raw_epsilon= epsilon; G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep(); G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();; if( epsilon < epsilonMin ) epsilon = epsilonMin; if( epsilon > epsilonMax ) epsilon = epsilonMax; SetEpsilonStep( epsilon ); // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon // << " final= " << epsilon << G4endl; // Shorten the proposed step in case of earlier problems (zero steps) // if( fNoZeroStep > fActionThreshold_NoZeroSteps ) { G4double stepTrial; stepTrial= fFull_CurveLen_of_LastAttempt; if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) ) stepTrial= fLast_ProposedStepLength; G4double decreaseFactor = 0.9; // Unused default if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps) && (stepTrial > 1000.0*kCarTolerance) ) { // Ensure quicker convergence // decreaseFactor= 0.1; } else { // We are in significant difficulties, probably at a boundary that // is either geometrically sharp or between very different materials. // Careful decreases to cope with tolerance are required. // if( stepTrial > 1000.0*kCarTolerance ) decreaseFactor = 0.25; // Try slow decreases else if( stepTrial > 100.0*kCarTolerance ) decreaseFactor= 0.5; // Try slower decreases else if( stepTrial > 10.0*kCarTolerance ) decreaseFactor= 0.75; // Try even slower decreases else decreaseFactor= 0.9; // Try very slow decreases } stepTrial *= decreaseFactor; #ifdef G4DEBUG_FIELD PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor, stepTrial, pFieldTrack); #endif if( stepTrial == 0.0 ) { G4cout << " G4PropagatorInField::ComputeStep " << " Particle abandoned due to lack of progress in field." << G4endl << " Properties : " << pFieldTrack << " " << G4endl; G4cerr << " G4PropagatorInField::ComputeStep " << " ERROR : attempting a zero step= " << stepTrial << G4endl << " while attempting to progress after " << fNoZeroStep << " trial steps. Will abandon step." << G4endl; fParticleIsLooping= true; return 0; // = stepTrial; } if( stepTrial < CurrentProposedStepLength ) CurrentProposedStepLength = stepTrial; } fLast_ProposedStepLength = CurrentProposedStepLength; G4int do_loop_count = 0; do { G4FieldTrack SubStepStartState = CurrentState; G4ThreeVector SubStartPoint = CurrentState.GetPosition(); if( !first_substep) { fNavigator->LocateGlobalPointWithinVolume( SubStartPoint ); } // How far to attempt to move the particle ! // h_TrialStepSize = CurrentProposedStepLength - StepTaken; // Integrate as far as "chord miss" rule allows. // s_length_taken = GetChordFinder()->AdvanceChordLimited( CurrentState, // Position & velocity h_TrialStepSize, fEpsilonStep, fPreviousSftOrigin, fPreviousSafety ); // CurrentState is now updated with the final position and velocity. fFull_CurveLen_of_LastAttempt = s_length_taken; G4ThreeVector EndPointB = CurrentState.GetPosition(); G4ThreeVector InterSectionPointE; G4double LinearStepLength; // Intersect chord AB with geometry intersects= IntersectChord( SubStartPoint, EndPointB, NewSafety, LinearStepLength, InterSectionPointE ); // E <- Intersection Point of chord AB and either volume A's surface // or a daughter volume's surface .. if( first_substep ) { currentSafety = NewSafety; } // Updating safety in other steps is potential future extention if( intersects ) { G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct // Find the intersection point of AB true path with the surface // of vol(A), if it exists. Start with point E as first "estimate". G4bool recalculatedEndPt= false; G4bool found_intersection = fIntersectionLocator-> EstimateIntersectionPoint( SubStepStartState, CurrentState, InterSectionPointE, IntersectPointVelct_G, recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin); intersects = intersects && found_intersection; if( found_intersection ) { End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ... StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength() - OriginalState.GetCurveLength(); } else { // intersects= false; // "Minor" chords do not intersect if( recalculatedEndPt ){ CurrentState= IntersectPointVelct_G; } } } if( !intersects ) { StepTaken += s_length_taken; // For smooth trajectory display (jacek 01/11/2002) if (fpTrajectoryFilter) { fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition()); } } first_substep = false; #ifdef G4DEBUG_FIELD if( fNoZeroStep > fActionThreshold_NoZeroSteps ) { printStatus( SubStepStartState, // or OriginalState, CurrentState, CurrentProposedStepLength, NewSafety, do_loop_count, pPhysVol ); } #endif #ifdef G4VERBOSE if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) { if( do_loop_count == fMax_loop_count-9 ){ G4cout << "G4PropagatorInField::ComputeStep " << " Difficult track - taking many sub steps." << G4endl; } printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength, NewSafety, do_loop_count, pPhysVol ); } #endif do_loop_count++; } while( (!intersects ) && (StepTaken + kCarTolerance < CurrentProposedStepLength) && ( do_loop_count < fMax_loop_count ) ); if( do_loop_count >= fMax_loop_count ) { fParticleIsLooping = true; if ( fVerboseLevel > 0 ){ G4cout << "G4PropagateInField: Killing looping particle " // << " of " << energy << " energy " << " after " << do_loop_count << " field substeps " << " totaling " << StepTaken / mm << " mm " ; if( pPhysVol ) G4cout << " in the volume " << pPhysVol->GetName() ; else G4cout << " in unknown or null volume. " ; G4cout << G4endl; } } if( !intersects ) { // Chord AB or "minor chords" do not intersect // B is the endpoint Step of the current Step. // End_PointAndTangent = CurrentState; TruePathLength = StepTaken; } // Set pFieldTrack to the return value // pFieldTrack = End_PointAndTangent; #ifdef G4VERBOSE // Check that "s" is correct // if( std::fabs(OriginalState.GetCurveLength() + TruePathLength - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength ) { G4cerr << " ERROR - G4PropagatorInField::ComputeStep():" << G4endl << " Curve length mis-match, is advancement wrong ? " << G4endl; G4cerr << " The curve length of the endpoint should be: " << OriginalState.GetCurveLength() + TruePathLength << G4endl << " and it is instead: " << End_PointAndTangent.GetCurveLength() << "." << G4endl << " A difference of: " << OriginalState.GetCurveLength() + TruePathLength - End_PointAndTangent.GetCurveLength() << G4endl; G4cerr << " Original state= " << OriginalState << G4endl << " Proposed state= " << End_PointAndTangent << G4endl; G4Exception("G4PropagatorInField::ComputeStep()", "IncorrectProposedEndPoint", FatalException, "Curve length mis-match between original state and proposed endpoint of propagation."); } #endif // In particular anomalous cases, we can get repeated zero steps // In order to correct this efficiently, we identify these cases // and only take corrective action when they occur. // if( TruePathLength < 0.5*kCarTolerance ) fNoZeroStep++; else fNoZeroStep = 0; if( fNoZeroStep > fAbandonThreshold_NoZeroSteps ) { fParticleIsLooping = true; G4cout << " WARNING - G4PropagatorInField::ComputeStep():" << G4endl << " Zero progress for " << fNoZeroStep << " attempted steps." << G4endl; G4cout << "Proposed Step is "<GetName() ; else G4cout << " in unknown or null volume. " ; G4cout << G4endl; if ( fVerboseLevel > 2 ) G4cout << " Particle that is stuck will be killed." << G4endl; fNoZeroStep = 0; } return TruePathLength; } /////////////////////////////////////////////////////////////////////////// // // Dumps status of propagator. void G4PropagatorInField::printStatus( const G4FieldTrack& StartFT, const G4FieldTrack& CurrentFT, G4double requestStep, G4double safety, G4int stepNo, G4VPhysicalVolume* startVolume) { const G4int verboseLevel= fVerboseLevel; const G4ThreeVector StartPosition = StartFT.GetPosition(); const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir(); const G4ThreeVector CurrentPosition = CurrentFT.GetPosition(); const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir(); G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength(); if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) ) { static G4int noPrecision= 4; G4cout.precision(noPrecision); // G4cout.setf(ios_base::fixed,ios_base::floatfield); G4cout << std::setw( 6) << " " << std::setw( 25) << " Current Position and Direction" << " " << G4endl; G4cout << std::setw( 5) << "Step#" << std::setw(10) << " s " << " " << std::setw(10) << "X(mm)" << " " << std::setw(10) << "Y(mm)" << " " << std::setw(10) << "Z(mm)" << " " << std::setw( 7) << " N_x " << " " << std::setw( 7) << " N_y " << " " << std::setw( 7) << " N_z " << " " ; // << G4endl; G4cout // << " >>> " << std::setw( 7) << " Delta|N|" << " " // << std::setw( 7) << " Delta(N_z) " << " " << std::setw( 9) << "StepLen" << " " << std::setw(12) << "StartSafety" << " " << std::setw( 9) << "PhsStep" << " "; if( startVolume ) { G4cout << std::setw(18) << "NextVolume" << " "; } G4cout << G4endl; } if((stepNo == 0) && (verboseLevel <=3)){ // Recurse to print the start values // printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume); } if( verboseLevel <= 3 ) { if( stepNo >= 0) G4cout << std::setw( 4) << stepNo << " "; else G4cout << std::setw( 5) << "Start" ; G4cout.precision(8); G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " "; G4cout.precision(8); G4cout << std::setw(10) << CurrentPosition.x() << " " << std::setw(10) << CurrentPosition.y() << " " << std::setw(10) << CurrentPosition.z() << " "; G4cout.precision(4); G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " " << std::setw( 7) << CurrentUnitVelocity.y() << " " << std::setw( 7) << CurrentUnitVelocity.z() << " "; // G4cout << G4endl; // G4cout << " >>> " ; G4cout.precision(3); G4cout << std::setw( 7) << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag() << " "; // << std::setw( 7) << CurrentUnitVelocity.z() - InitialUnitVelocity.z() << " "; G4cout << std::setw( 9) << step_len << " "; G4cout << std::setw(12) << safety << " "; if( requestStep != -1.0 ) G4cout << std::setw( 9) << requestStep << " "; else G4cout << std::setw( 9) << "Init/NotKnown" << " "; if( startVolume != 0) { G4cout << std::setw(12) << startVolume->GetName() << " "; } #if 0 else { if( step_len != -1 ) G4cout << std::setw(12) << "OutOfWorld" << " "; else G4cout << std::setw(12) << "NotGiven" << " "; } #endif G4cout << G4endl; } else // if( verboseLevel > 3 ) { // Multi-line output G4cout << "Step taken was " << step_len << " out of PhysicalStep= " << requestStep << G4endl; G4cout << "Final safety is: " << safety << G4endl; G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag() << G4endl; G4cout << G4endl; } } /////////////////////////////////////////////////////////////////////////// // // Prints Step diagnostics void G4PropagatorInField::PrintStepLengthDiagnostic( G4double CurrentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack& ) { G4cout << " PiF: NoZeroStep= " << fNoZeroStep << " CurrentProposedStepLength= " << CurrentProposedStepLength << " Full_curvelen_last=" << fFull_CurveLen_of_LastAttempt << " last proposed step-length= " << fLast_ProposedStepLength << " decreate factor = " << decreaseFactor << " step trial = " << stepTrial << G4endl; } // Access the points which have passed through the filter. The // points are stored as ThreeVectors for the initial impelmentation // only (jacek 30/10/2002) // Responsibility for deleting the points lies with // SmoothTrajectoryPoint, which is the points' final // destination. The points pointer is set to NULL, to ensure that // the points are not re-used in subsequent steps, therefore THIS // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002) std::vector* G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const { // NB, GimmeThePointsAndForgetThem really forgets them, so it can // only be called (exactly) once for each step. if (fpTrajectoryFilter) { return fpTrajectoryFilter->GimmeThePointsAndForgetThem(); } else { return 0; } } void G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter) { fpTrajectoryFilter = filter; } void G4PropagatorInField::ClearPropagatorState() { // Goal: Clear all memory of previous steps, cached information fParticleIsLooping= false; fNoZeroStep= 0; End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.), G4ThreeVector(0.,0.,0.), 0.0,0.0,0.0,0.0,0.0); fFull_CurveLen_of_LastAttempt = -1; fLast_ProposedStepLength = -1; fPreviousSftOrigin= G4ThreeVector(0.,0.,0.); fPreviousSafety= 0.0; } G4FieldManager* G4PropagatorInField:: FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume) { G4FieldManager* currentFieldMgr; currentFieldMgr = fDetectorFieldMgr; if( pCurrentPhysicalVolume) { G4FieldManager *newFieldMgr = 0; newFieldMgr= pCurrentPhysicalVolume->GetLogicalVolume()->GetFieldManager(); if ( newFieldMgr ) currentFieldMgr = newFieldMgr; } fCurrentFieldMgr= currentFieldMgr; // Flag that field manager has been set. fSetFieldMgr= true; return currentFieldMgr; } G4int G4PropagatorInField::SetVerboseLevel( G4int level ) { G4int oldval= fVerboseLevel; fVerboseLevel= level; // Forward the verbose level 'reduced' to ChordFinder, // MagIntegratorDriver ... ? // G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver(); integrDriver->SetVerboseLevel( fVerboseLevel - 2 ); G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl; return oldval; }