G4ChordFinderSaf.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 #include "G4ChordFinderSaf.hh"
00029 #include <iomanip>
00030 
00031 // ..........................................................................
00032 
00033 G4ChordFinderSaf::G4ChordFinderSaf(G4MagInt_Driver* pIntegrationDriver)
00034   : G4ChordFinder(pIntegrationDriver)
00035 {
00036     // check the values and set the other parameters
00037     // fNoInitialRadBig=0;    fNoInitialRadSmall=0;  
00038     // fNoTrialsRadBig=0;     fNoTrialsRadSmall=0; 
00039    
00040 }
00041 
00042 // ..........................................................................
00043 
00044 G4ChordFinderSaf::G4ChordFinderSaf( G4MagneticField*        theMagField,
00045                               G4double                stepMinimum, 
00046                               G4MagIntegratorStepper* pItsStepper )
00047   : G4ChordFinder( theMagField, stepMinimum, pItsStepper )
00048 {
00049   //  Let  G4ChordFinder create the  Driver, the Stepper and EqRhs ...
00050   // ...
00051 }
00052 
00053 // ......................................................................
00054 
00055 // ......................................................................
00056 
00057 G4ChordFinderSaf::~G4ChordFinderSaf()
00058 {
00059   if( SetVerbose(0) ) { PrintStatistics();  } 
00060    // Set verbosity 0, so that will be called in base class again
00061 }
00062 
00063 void
00064 G4ChordFinderSaf::PrintStatistics()
00065 {
00066   // Print Statistics
00067   G4cout << "G4ChordFinderSaf statistics report: " << G4endl;
00068   G4ChordFinder::PrintStatistics();
00069 
00070 /*******************
00071   G4cout 
00072     << "  No radbig calls " << std::setw(10) << fNoInitialRadBig 
00073     << " trials " << std::setw(10) << fNoTrialsRadBig
00074     << "  - over " << std::setw(10) << fNoTrialsRadBig - fNoInitialRadBig
00075     << G4endl
00076     << "  No radsm  calls " << std::setw(10) << fNoInitialRadSmall 
00077     << " trials "  << std::setw(10) << fNoTrialsRadSmall
00078     << "  - over " << std::setw(10) << fNoTrialsRadSmall - fNoInitialRadSmall
00079     << G4endl;
00080   G4cout
00081     << "  *** Limiting stepTrial via if Delta_chord < R_curvature " 
00082     << "   for large to angle from Delta_chord / R_curv "  
00083     << "   and for small with multiple " << GetMultipleRadius()
00084     << G4endl; 
00085 ********************/
00086 }
00087 
00088 
00089 // G4SafetyDist::
00090 // inline
00091 G4double
00092 CalculatePointSafety(G4ThreeVector safetyOrigin,
00093                      G4double      safetyRadius,
00094                      G4ThreeVector point)
00095 {
00096   G4double      pointSafety= 0.0; 
00097 
00098   G4ThreeVector OriginShift = point - safetyOrigin ;
00099   G4double      MagSqShift  = OriginShift.mag2() ;
00100   if( MagSqShift < sqr(safetyRadius) ){ 
00101     pointSafety = safetyRadius - std::sqrt(MagSqShift) ;  
00102   }
00103 
00104   return pointSafety;
00105 }
00106 
00107 // inline
00108 G4bool
00109 CalculatePointInside(G4ThreeVector safetyOrigin,
00110                      G4double      safetyRadius,
00111                      G4ThreeVector point)
00112 {
00113   G4ThreeVector OriginShift = point - safetyOrigin ;
00114   return ( OriginShift.mag2() < safetyRadius*safetyRadius ); 
00115 }
00116 
00117 G4double
00118 G4ChordFinderSaf::FindNextChord( const  G4FieldTrack&  yStart,
00119                                      G4double     stepMax,
00120                                      G4FieldTrack&   yEnd, // Endpoint
00121                                      G4double&   dyErrPos, // Error of endpoint
00122                                      G4double    epsStep,
00123                                      G4double*  pStepForAccuracy, 
00124                                     const G4ThreeVector latestSafetyOrigin,
00125                                     G4double       latestSafetyRadius 
00126                                         )
00127   // Returns Length of Step taken
00128 {
00129   // G4int       stepRKnumber=0;
00130   G4FieldTrack yCurrent=  yStart;  
00131   G4double    stepTrial, stepForAccuracy;
00132   G4double    dydx[G4FieldTrack::ncompSVEC]; 
00133 
00134   //  1.)  Try to "leap" to end of interval
00135   //  2.)  Evaluate if resulting chord gives d_chord that is good enough.
00136   //     2a.)  If d_chord is not good enough, find one that is.
00137   
00138   G4bool    validEndPoint= false;
00139   G4double  dChordStep, lastStepLength; //  stepOfLastGoodChord;
00140 
00141   GetIntegrationDriver()-> GetDerivatives( yCurrent, dydx )  ;
00142 
00143   G4int     noTrials=0;
00144   const G4double safetyFactor= GetFirstFraction(); // was 0.999
00145 
00146   // Figure out the starting safety
00147   G4double      startSafety=  
00148     CalculatePointSafety( latestSafetyOrigin, 
00149                           latestSafetyRadius, 
00150                           yCurrent.GetPosition() );
00151 
00152   G4double
00153     likelyGood = std::max( startSafety , 
00154                         safetyFactor *  GetLastStepEstimateUnc() );
00155 
00156   stepTrial  = std::min( stepMax,  likelyGood ); 
00157 
00158   G4MagInt_Driver *pIntgrDriver= G4ChordFinder::GetIntegrationDriver(); 
00159   G4double newStepEst_Uncons= 0.0; 
00160   do
00161   { 
00162      G4double stepForChord;  
00163      yCurrent = yStart;    // Always start from initial point
00164 
00165    //            ************
00166      pIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial, 
00167                                  dChordStep, dyErrPos);
00168      //            ************
00169 
00170      G4ThreeVector EndPointCand= yCurrent.GetPosition(); 
00171      G4bool  endPointInSafetySphere= 
00172        CalculatePointInside(latestSafetyOrigin, latestSafetyRadius, EndPointCand);
00173 
00174      // We check whether the criterion is met here.
00175      validEndPoint = AcceptableMissDist(dChordStep)
00176                        || endPointInSafetySphere;
00177                       //  && (dyErrPos < eps) ;
00178 
00179      lastStepLength = stepTrial; 
00180 
00181      // This method estimates to step size for a good chord.
00182      stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
00183 
00184      if( ! validEndPoint ) {
00185         if( stepTrial<=0.0 )
00186           stepTrial = stepForChord; 
00187         else if (stepForChord <= stepTrial) 
00188           // Reduce by a fraction, possibly up to 20% 
00189           stepTrial = std::min( stepForChord, 
00190                                 GetFractionLast() * stepTrial); 
00191         else
00192           stepTrial *= 0.1;
00193 
00194         // if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
00195      }
00196      // #ifdef  TEST_CHORD_PRINT
00197      // TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
00198      // #endif
00199 
00200      noTrials++; 
00201   }
00202   while( ! validEndPoint );   // End of do-while  RKD 
00203 
00204   AccumulateStatistics( noTrials );
00205 
00206   // Should we update newStepEst_Uncons for a 'long step' via safety ??
00207   if( newStepEst_Uncons > 0.0  ){ 
00208     SetLastStepEstimateUnc( newStepEst_Uncons );
00209   }
00210 
00211   // stepOfLastGoodChord = stepTrial;
00212 
00213   if( pStepForAccuracy ){ 
00214      // Calculate the step size required for accuracy, if it is needed
00215      G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
00216      if( dyErr_relative > 1.0 ) {
00217         stepForAccuracy =
00218            pIntgrDriver->ComputeNewStepSize( dyErr_relative,
00219                                              lastStepLength );
00220      }else{
00221         stepForAccuracy = 0.0;   // Convention to show step was ok 
00222      }
00223      *pStepForAccuracy = stepForAccuracy;
00224   }
00225 
00226 #ifdef  TEST_CHORD_PRINT
00227   static int dbg=0;
00228   if( dbg ) 
00229     G4cout << "ChordF/FindNextChord:  NoTrials= " << noTrials 
00230            << " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
00231 #endif
00232 
00233   yEnd=  yCurrent;  
00234   return stepTrial; 
00235 }

Generated on Mon May 27 17:47:53 2013 for Geant4 by  doxygen 1.4.7