G4PropagatorInField.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // 
00028 // 
00029 //  This class implements an algorithm to track a particle in a
00030 //  non-uniform magnetic field. It utilises an ODE solver (with
00031 //  the Runge - Kutta method) to evolve the particle, and drives it
00032 //  until the particle has traveled a set distance or it enters a new 
00033 //  volume.
00034 //                                                                     
00035 // 14.10.96 John Apostolakis,   design and implementation
00036 // 17.03.97 John Apostolakis,   renaming new set functions being added
00037 //
00038 // $Id$
00039 // GEANT4 tag $ Name:  $
00040 // ---------------------------------------------------------------------------
00041 
00042 #include <iomanip>
00043 
00044 #include "G4PropagatorInField.hh"
00045 #include "G4ios.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4ThreeVector.hh"
00048 #include "G4VPhysicalVolume.hh"
00049 #include "G4Navigator.hh"
00050 #include "G4GeometryTolerance.hh"
00051 #include "G4VCurvedTrajectoryFilter.hh"
00052 #include "G4ChordFinder.hh"
00053 #include "G4MultiLevelLocator.hh"
00054 
00056 //
00057 // Constructors and destructor
00058 
00059 G4PropagatorInField::G4PropagatorInField( G4Navigator    *theNavigator, 
00060                                           G4FieldManager *detectorFieldMgr,
00061                                           G4VIntersectionLocator *vLocator  )
00062   : 
00063     fMax_loop_count(1000),
00064     fUseSafetyForOptimisation(true),   // (false) is less sensitive to incorrect safety
00065     fZeroStepThreshold( 0.0 ),         // length of what is recognised as 'zero' step
00066     fDetectorFieldMgr(detectorFieldMgr), 
00067     fpTrajectoryFilter( 0 ),
00068     fNavigator(theNavigator),
00069     fCurrentFieldMgr(detectorFieldMgr),
00070     fSetFieldMgr(false),
00071     fCharge(0.0), fInitialMomentumModulus(0.0), fMass(0.0),
00072     End_PointAndTangent(G4ThreeVector(0.,0.,0.),
00073                         G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0),
00074     fParticleIsLooping(false),
00075     fNoZeroStep(0), 
00076     fVerboseLevel(0)
00077 {
00078   if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
00079   else                  { fEpsilonStep= 1.0e-5; } 
00080   fActionThreshold_NoZeroSteps = 2; 
00081   fSevereActionThreshold_NoZeroSteps = 10; 
00082   fAbandonThreshold_NoZeroSteps = 50; 
00083   fFull_CurveLen_of_LastAttempt = -1; 
00084   fLast_ProposedStepLength = -1;
00085   fLargestAcceptableStep = 1000.0 * meter;
00086 
00087   fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
00088   fPreviousSafety= 0.0;
00089   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00090   fZeroStepThreshold= std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer );
00091 
00092 #ifdef G4DEBUG_FIELD
00093   G4cout << " PiF: Zero Step Threshold set to "
00094          << fZeroStepThreshold / millimeter
00095          << " mm." << G4endl;
00096   G4cout << " PiF:   Value of kCarTolerance = "
00097          << kCarTolerance / millimeter 
00098          << " mm. " << G4endl;
00099 #endif 
00100 
00101   // Definding Intersection Locator and his parameters
00102   if(vLocator==0){
00103     fIntersectionLocator= new G4MultiLevelLocator(theNavigator);
00104     fAllocatedLocator=true;
00105   }else{
00106     fIntersectionLocator=vLocator;
00107     fAllocatedLocator=false;
00108   }
00109   RefreshIntersectionLocator();  //  Copy all relevant parameters 
00110 }
00111 
00112 G4PropagatorInField::~G4PropagatorInField()
00113 {
00114   if(fAllocatedLocator)delete  fIntersectionLocator; 
00115 }
00116 
00117 // Update the IntersectionLocator with current parameters
00118 void
00119 G4PropagatorInField::RefreshIntersectionLocator()
00120 {
00121   fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
00122   fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
00123   fIntersectionLocator->SetChordFinderFor(GetChordFinder());
00124   fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
00125 }
00127 //
00128 // Compute the next geometric Step
00129 
00130 G4double
00131 G4PropagatorInField::ComputeStep(
00132                 G4FieldTrack&      pFieldTrack,
00133                 G4double           CurrentProposedStepLength,
00134                 G4double&          currentSafety,                // IN/OUT
00135                 G4VPhysicalVolume* pPhysVol)
00136 {  
00137   // If CurrentProposedStepLength is too small for finding Chords
00138   // then return with no action (for now - TODO: some action)
00139   //
00140   if(CurrentProposedStepLength<kCarTolerance)
00141   {
00142     return kInfinity;
00143   }
00144 
00145   // Introducing smooth trajectory display (jacek 01/11/2002)
00146   //
00147   if (fpTrajectoryFilter)
00148   {
00149     fpTrajectoryFilter->CreateNewTrajectorySegment();
00150   }
00151 
00152   // Parameters for adaptive Runge-Kutta integration
00153   
00154   G4double      h_TrialStepSize;        // 1st Step Size 
00155   G4double      TruePathLength = CurrentProposedStepLength;
00156   G4double      StepTaken = 0.0; 
00157   G4double      s_length_taken, epsilon ; 
00158   G4bool        intersects;
00159   G4bool        first_substep = true;
00160 
00161   G4double      NewSafety;
00162   fParticleIsLooping = false;
00163 
00164   // If not yet done, 
00165   //   Set the field manager to the local  one if the volume has one, 
00166   //                      or to the global one if not
00167   //
00168   if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol ); 
00169   // For the next call, the field manager must again be set
00170   fSetFieldMgr= false;
00171 
00172   GetChordFinder()->SetChargeMomentumMass(fCharge,fInitialMomentumModulus,fMass);
00173 
00174  // Values for Intersection Locator has to be updated on each call for the
00175  // case that CurrentFieldManager has changed from the one of previous step
00176   RefreshIntersectionLocator();
00177 
00178   G4FieldTrack  CurrentState(pFieldTrack);
00179   G4FieldTrack  OriginalState = CurrentState;
00180 
00181   // If the Step length is "infinite", then an approximate-maximum Step
00182   // length (used to calculate the relative accuracy) must be guessed.
00183   //
00184   if( CurrentProposedStepLength >= fLargestAcceptableStep )
00185   {
00186     G4ThreeVector StartPointA, VelocityUnit;
00187     StartPointA  = pFieldTrack.GetPosition();
00188     VelocityUnit = pFieldTrack.GetMomentumDir();
00189 
00190     G4double trialProposedStep = 1.e2 * ( 10.0 * cm + 
00191       fNavigator->GetWorldVolume()->GetLogicalVolume()->
00192                   GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
00193     CurrentProposedStepLength= std::min( trialProposedStep,
00194                                            fLargestAcceptableStep ); 
00195   }
00196   epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
00197   // G4double raw_epsilon= epsilon;
00198   G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
00199   G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();; 
00200   if( epsilon < epsilonMin ) epsilon = epsilonMin;
00201   if( epsilon > epsilonMax ) epsilon = epsilonMax;
00202   SetEpsilonStep( epsilon );
00203 
00204   // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
00205   //        << " final= " << epsilon << G4endl;
00206 
00207   //  Shorten the proposed step in case of earlier problems (zero steps)
00208   // 
00209   if( fNoZeroStep > fActionThreshold_NoZeroSteps )
00210   {
00211     G4double stepTrial;
00212 
00213     stepTrial= fFull_CurveLen_of_LastAttempt; 
00214     if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) ) 
00215       stepTrial= fLast_ProposedStepLength; 
00216 
00217     G4double decreaseFactor = 0.9; // Unused default
00218     if(   (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
00219        && (stepTrial > 100.0*fZeroStepThreshold) )
00220     {
00221       // Attempt quick convergence
00222       //
00223       decreaseFactor= 0.25;
00224     } 
00225     else
00226     {
00227       // We are in significant difficulties, probably at a boundary that
00228       // is either geometrically sharp or between very different materials.
00229       // Careful decreases to cope with tolerance are required.
00230       //
00231       if( stepTrial > 100.0*fZeroStepThreshold )
00232         decreaseFactor = 0.35;     // Try decreasing slower
00233       else if( stepTrial > 30.0*fZeroStepThreshold )
00234         decreaseFactor= 0.5;       // Try yet slower decrease
00235       else if( stepTrial > 10.0*fZeroStepThreshold )
00236         decreaseFactor= 0.75;      // Try even slower decreases
00237       else
00238         decreaseFactor= 0.9;       // Try very slow decreases
00239      }
00240      stepTrial *= decreaseFactor;
00241 
00242 #ifdef G4DEBUG_FIELD
00243      G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
00244             << "  Decreasing step -  " 
00245             << " decreaseFactor= " << std::setw(8) << decreaseFactor 
00246             << " stepTrial = "     << std::setw(18) << stepTrial << " "
00247             << " fZeroStepThreshold = " << fZeroStepThreshold << G4endl;
00248      PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
00249                                stepTrial, pFieldTrack);
00250 #endif
00251      if( stepTrial == 0.0 )  //  Change to make it < 0.1 * kCarTolerance ??
00252      {
00253        std::ostringstream message;
00254        message << "Particle abandoned due to lack of progress in field."
00255                << G4endl
00256                << "  Properties : " << pFieldTrack << G4endl
00257                << "  Attempting a zero step = " << stepTrial << G4endl
00258                << "  while attempting to progress after " << fNoZeroStep
00259                << " trial steps. Will abandon step.";
00260        G4Exception("G4PropagatorInField::ComputeStep()", "GeomNav1002",
00261                    JustWarning, message);
00262        fParticleIsLooping= true;
00263        return 0;  // = stepTrial;
00264      }
00265      if( stepTrial < CurrentProposedStepLength )
00266        CurrentProposedStepLength = stepTrial;
00267   }
00268   fLast_ProposedStepLength = CurrentProposedStepLength;
00269 
00270   G4int do_loop_count = 0; 
00271   do
00272   { 
00273     G4FieldTrack SubStepStartState = CurrentState;
00274     G4ThreeVector SubStartPoint = CurrentState.GetPosition(); 
00275 
00276     if( !first_substep) {
00277       fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
00278     }
00279 
00280     // How far to attempt to move the particle !
00281     //
00282     h_TrialStepSize = CurrentProposedStepLength - StepTaken;
00283 
00284     // Integrate as far as "chord miss" rule allows.
00285     //
00286     s_length_taken = GetChordFinder()->AdvanceChordLimited( 
00287                              CurrentState,    // Position & velocity
00288                              h_TrialStepSize,
00289                              fEpsilonStep,
00290                              fPreviousSftOrigin,
00291                              fPreviousSafety
00292                              );
00293     //  CurrentState is now updated with the final position and velocity. 
00294 
00295     fFull_CurveLen_of_LastAttempt = s_length_taken;
00296 
00297     G4ThreeVector  EndPointB = CurrentState.GetPosition(); 
00298     G4ThreeVector  InterSectionPointE;
00299     G4double       LinearStepLength;
00300  
00301     // Intersect chord AB with geometry
00302     intersects= IntersectChord( SubStartPoint, EndPointB, 
00303                                 NewSafety,     LinearStepLength, 
00304                                 InterSectionPointE );
00305     // E <- Intersection Point of chord AB and either volume A's surface 
00306     //                                  or a daughter volume's surface ..
00307 
00308     if( first_substep ) { 
00309        currentSafety = NewSafety;
00310     } // Updating safety in other steps is potential future extention
00311 
00312     if( intersects )
00313     {
00314        G4FieldTrack IntersectPointVelct_G(CurrentState);  // FT-Def-Construct
00315 
00316        // Find the intersection point of AB true path with the surface
00317        //   of vol(A), if it exists. Start with point E as first "estimate".
00318        G4bool recalculatedEndPt= false;
00319        
00320          G4bool found_intersection = fIntersectionLocator->
00321          EstimateIntersectionPoint( SubStepStartState, CurrentState, 
00322                                   InterSectionPointE, IntersectPointVelct_G,
00323                                   recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
00324        intersects = intersects && found_intersection;
00325        if( found_intersection ) {        
00326           End_PointAndTangent= IntersectPointVelct_G;  // G is our EndPoint ...
00327           StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
00328                                       - OriginalState.GetCurveLength();
00329        } else {
00330           // intersects= false;          // "Minor" chords do not intersect
00331           if( recalculatedEndPt ){
00332              CurrentState= IntersectPointVelct_G; 
00333           }
00334        }
00335     }
00336     if( !intersects )
00337     {
00338       StepTaken += s_length_taken; 
00339       // For smooth trajectory display (jacek 01/11/2002)
00340       if (fpTrajectoryFilter) {
00341         fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
00342       }
00343     }
00344     first_substep = false;
00345 
00346 #ifdef G4DEBUG_FIELD
00347     if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
00348       printStatus( SubStepStartState,  // or OriginalState,
00349                    CurrentState,  CurrentProposedStepLength, 
00350                    NewSafety,     do_loop_count,  pPhysVol );
00351     }
00352     if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
00353       if( do_loop_count == fMax_loop_count-9 ){
00354         G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
00355                << "  Difficult track - taking many sub steps." << G4endl;
00356       }
00357       printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength, 
00358                    NewSafety, do_loop_count, pPhysVol );
00359     }
00360 #endif
00361 
00362     do_loop_count++;
00363 
00364   } while( (!intersects )
00365         && (StepTaken + kCarTolerance < CurrentProposedStepLength)  
00366         && ( do_loop_count < fMax_loop_count ) );
00367 
00368   if( do_loop_count >= fMax_loop_count  )
00369   {
00370     fParticleIsLooping = true;
00371 
00372     if ( fVerboseLevel > 0 )
00373     {
00374        G4cout << " G4PropagateInField::ComputeStep(): " << G4endl
00375               << "  Killing looping particle " 
00376               // << " of " << energy  << " energy "
00377               << " after " << do_loop_count << " field substeps "
00378               << " totaling " << StepTaken / mm << " mm " ;
00379        if( pPhysVol )
00380           G4cout << " in volume " << pPhysVol->GetName() ; 
00381        else
00382          G4cout << " in unknown or null volume. " ; 
00383        G4cout << G4endl;
00384     }
00385   }
00386 
00387   if( !intersects )
00388   {
00389     // Chord AB or "minor chords" do not intersect
00390     // B is the endpoint Step of the current Step.
00391     //
00392     End_PointAndTangent = CurrentState; 
00393     TruePathLength = StepTaken;
00394   }
00395   
00396   // Set pFieldTrack to the return value
00397   //
00398   pFieldTrack = End_PointAndTangent;
00399 
00400 #ifdef G4VERBOSE
00401   // Check that "s" is correct
00402   //
00403   if( std::fabs(OriginalState.GetCurveLength() + TruePathLength 
00404       - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
00405   {
00406     std::ostringstream message;
00407     message << "Curve length mis-match between original state "
00408             << "and proposed endpoint of propagation." << G4endl
00409             << "  The curve length of the endpoint should be: " 
00410             << OriginalState.GetCurveLength() + TruePathLength << G4endl
00411             << "  and it is instead: "
00412             << End_PointAndTangent.GetCurveLength() << "." << G4endl
00413             << "  A difference of: "
00414             << OriginalState.GetCurveLength() + TruePathLength 
00415                - End_PointAndTangent.GetCurveLength() << G4endl
00416             << "  Original state = " << OriginalState   << G4endl
00417             << "  Proposed state = " << End_PointAndTangent;
00418     G4Exception("G4PropagatorInField::ComputeStep()",
00419                 "GeomNav0003", FatalException, message);
00420   }
00421 #endif
00422 
00423   // In particular anomalous cases, we can get repeated zero steps
00424   // In order to correct this efficiently, we identify these cases
00425   // and only take corrective action when they occur.
00426   // 
00427   if( ( (TruePathLength < fZeroStepThreshold) 
00428         && ( TruePathLength+kCarTolerance < CurrentProposedStepLength  ) 
00429         ) 
00430       || ( TruePathLength < 0.5*kCarTolerance )
00431     )
00432   {
00433     fNoZeroStep++;
00434   }
00435   else{
00436     fNoZeroStep = 0;
00437   }
00438 
00439   if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
00440   { 
00441      fParticleIsLooping = true;
00442      std::ostringstream message;
00443      message << "Particle is stuck; it will be killed." << G4endl
00444              << "  Zero progress for "  << fNoZeroStep << " attempted steps." 
00445              << G4endl
00446              << "  Proposed Step is " << CurrentProposedStepLength
00447              << " but Step Taken is "<< fFull_CurveLen_of_LastAttempt << G4endl
00448              << "  For Particle with Charge = " << fCharge
00449              << " Momentum = "<< fInitialMomentumModulus
00450              << " Mass = " << fMass << G4endl;
00451      if( pPhysVol )
00452        message << " in volume " << pPhysVol->GetName() ; 
00453      else
00454        message << " in unknown or null volume. " ; 
00455      G4Exception("G4PropagatorInField::ComputeStep()",
00456                  "GeomNav1002", JustWarning, message);
00457      fNoZeroStep = 0; 
00458   }
00459  
00460   return TruePathLength;
00461 }
00462 
00464 //
00465 // Dumps status of propagator.
00466 
00467 void
00468 G4PropagatorInField::printStatus( const G4FieldTrack&        StartFT,
00469                                   const G4FieldTrack&        CurrentFT, 
00470                                         G4double             requestStep, 
00471                                         G4double             safety,
00472                                         G4int                stepNo, 
00473                                         G4VPhysicalVolume*   startVolume)
00474 {
00475   const G4int verboseLevel=fVerboseLevel;
00476   const G4ThreeVector StartPosition       = StartFT.GetPosition();
00477   const G4ThreeVector StartUnitVelocity   = StartFT.GetMomentumDir();
00478   const G4ThreeVector CurrentPosition     = CurrentFT.GetPosition();
00479   const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
00480 
00481   G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
00482 
00483   G4int oldprec;   // cout/cerr precision settings
00484       
00485   if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
00486   {
00487     oldprec = G4cout.precision(4);
00488     G4cout << std::setw( 6)  << " " 
00489            << std::setw( 25) << " Current Position  and  Direction" << " "
00490            << G4endl; 
00491     G4cout << std::setw( 5) << "Step#" 
00492            << std::setw(10) << "  s  " << " "
00493            << std::setw(10) << "X(mm)" << " "
00494            << std::setw(10) << "Y(mm)" << " "  
00495            << std::setw(10) << "Z(mm)" << " "
00496            << std::setw( 7) << " N_x " << " "
00497            << std::setw( 7) << " N_y " << " "
00498            << std::setw( 7) << " N_z " << " " ;
00499     G4cout << std::setw( 7) << " Delta|N|" << " "
00500            << std::setw( 9) << "StepLen" << " "  
00501            << std::setw(12) << "StartSafety" << " "  
00502            << std::setw( 9) << "PhsStep" << " ";  
00503     if( startVolume )
00504       { G4cout << std::setw(18) << "NextVolume" << " "; }
00505     G4cout.precision(oldprec);
00506     G4cout << G4endl;
00507   }
00508   if((stepNo == 0) && (verboseLevel <=3))
00509   {
00510     // Recurse to print the start values
00511     //
00512     printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
00513   }
00514   if( verboseLevel <= 3 )
00515   {
00516     if( stepNo >= 0)
00517       { G4cout << std::setw( 4) << stepNo << " "; }
00518     else
00519       { G4cout << std::setw( 5) << "Start" ; }
00520     oldprec = G4cout.precision(8);
00521     G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " "; 
00522     G4cout.precision(8);
00523     G4cout << std::setw(10) << CurrentPosition.x() << " "
00524            << std::setw(10) << CurrentPosition.y() << " "
00525            << std::setw(10) << CurrentPosition.z() << " ";
00526     G4cout.precision(4);
00527     G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
00528            << std::setw( 7) << CurrentUnitVelocity.y() << " "
00529            << std::setw( 7) << CurrentUnitVelocity.z() << " ";
00530     G4cout.precision(3); 
00531     G4cout << std::setw( 7)
00532            << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " "; 
00533     G4cout << std::setw( 9) << step_len << " "; 
00534     G4cout << std::setw(12) << safety << " ";
00535     if( requestStep != -1.0 ) 
00536       { G4cout << std::setw( 9) << requestStep << " "; }
00537     else
00538       { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
00539     if( startVolume != 0)
00540       { G4cout << std::setw(12) << startVolume->GetName() << " "; }
00541     G4cout.precision(oldprec);
00542     G4cout << G4endl;
00543   }
00544   else // if( verboseLevel > 3 )
00545   {
00546     //  Multi-line output
00547       
00548     G4cout << "Step taken was " << step_len  
00549            << " out of PhysicalStep = " <<  requestStep << G4endl;
00550     G4cout << "Final safety is: " << safety << G4endl;
00551     G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
00552            << G4endl;
00553     G4cout << G4endl; 
00554   }
00555 }
00556 
00558 //
00559 // Prints Step diagnostics
00560 
00561 void 
00562 G4PropagatorInField::PrintStepLengthDiagnostic(
00563                           G4double CurrentProposedStepLength,
00564                           G4double decreaseFactor,
00565                           G4double stepTrial,
00566                     const G4FieldTrack& )
00567 {
00568   G4int  iprec= G4cout.precision(8); 
00569   G4cout << " " << std::setw(12) << " PiF: NoZeroStep " 
00570          << " " << std::setw(20) << " CurrentProposed len " 
00571          << " " << std::setw(18) << " Full_curvelen_last" 
00572          << " " << std::setw(18) << " last proposed len " 
00573          << " " << std::setw(18) << " decrease factor   " 
00574          << " " << std::setw(15) << " step trial  " 
00575          << G4endl;
00576 
00577   G4cout << " " << std::setw(10) << fNoZeroStep << "  "
00578          << " " << std::setw(20) << CurrentProposedStepLength
00579          << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
00580          << " " << std::setw(18) << fLast_ProposedStepLength 
00581          << " " << std::setw(18) << decreaseFactor
00582          << " " << std::setw(15) << stepTrial
00583          << G4endl;
00584   G4cout.precision( iprec ); 
00585 
00586 }
00587 
00588 // Access the points which have passed through the filter. The
00589 // points are stored as ThreeVectors for the initial impelmentation
00590 // only (jacek 30/10/2002)
00591 // Responsibility for deleting the points lies with
00592 // SmoothTrajectoryPoint, which is the points' final
00593 // destination. The points pointer is set to NULL, to ensure that
00594 // the points are not re-used in subsequent steps, therefore THIS
00595 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
00596 
00597 std::vector<G4ThreeVector>*
00598 G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
00599 {
00600   // NB, GimmeThePointsAndForgetThem really forgets them, so it can
00601   // only be called (exactly) once for each step.
00602 
00603   if (fpTrajectoryFilter)
00604   {
00605     return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
00606   }
00607   else
00608   {
00609     return 0;
00610   }
00611 }
00612 
00613 void 
00614 G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
00615 {
00616   fpTrajectoryFilter = filter;
00617 }
00618 
00619 void G4PropagatorInField::ClearPropagatorState()
00620 {
00621   // Goal: Clear all memory of previous steps,  cached information
00622 
00623   fParticleIsLooping= false;
00624   fNoZeroStep= 0;
00625 
00626   End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
00627                                      G4ThreeVector(0.,0.,0.),
00628                                      0.0,0.0,0.0,0.0,0.0); 
00629   fFull_CurveLen_of_LastAttempt = -1; 
00630   fLast_ProposedStepLength = -1;
00631 
00632   fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
00633   fPreviousSafety= 0.0;
00634 }
00635 
00636 G4FieldManager* G4PropagatorInField::
00637 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
00638 {
00639   G4FieldManager* currentFieldMgr;
00640 
00641   currentFieldMgr = fDetectorFieldMgr;
00642   if( pCurrentPhysicalVolume)
00643   {
00644      G4FieldManager *pRegionFieldMgr= 0, *localFieldMgr = 0;
00645      G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
00646 
00647      if( pLogicalVol ) { 
00648         // Value for Region, if any, Overrides 
00649         G4Region*  pRegion= pLogicalVol->GetRegion();
00650         if( pRegion ) { 
00651            pRegionFieldMgr= pRegion->GetFieldManager();
00652            if( pRegionFieldMgr ) 
00653              currentFieldMgr= pRegionFieldMgr;
00654         }
00655 
00656         // 'Local' Value from logical volume, if any, Overrides 
00657         localFieldMgr= pLogicalVol->GetFieldManager();
00658         if ( localFieldMgr ) 
00659            currentFieldMgr = localFieldMgr;
00660      }
00661   }
00662   fCurrentFieldMgr= currentFieldMgr;
00663 
00664   // Flag that field manager has been set.
00665   fSetFieldMgr= true;
00666 
00667   return currentFieldMgr;
00668 }
00669 
00670 G4int G4PropagatorInField::SetVerboseLevel( G4int level )
00671 {
00672   G4int oldval= fVerboseLevel;
00673   fVerboseLevel= level;
00674 
00675   // Forward the verbose level 'reduced' to ChordFinder,
00676   // MagIntegratorDriver ... ? 
00677   //
00678   G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver(); 
00679   integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
00680   G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
00681 
00682   return oldval;
00683 }

Generated on Mon May 27 17:49:26 2013 for Geant4 by  doxygen 1.4.7