#include <TrackPropagation/RungeKutta/interface/RK4PreciseStep.h>
Public Member Functions | |
double | distance (const CartesianState &a, const CartesianState &b) const |
CartesianState | operator() (const CartesianState &start, const RKCartesianDerivative &deriv, double step, double eps) const |
std::pair< CartesianState, double > | stepWithAccuracy (const CartesianState &start, const RKCartesianDerivative &deriv, double step) const |
Private Member Functions | |
bool | verbose () const |
Definition at line 9 of file RK4PreciseStep.h.
double RK4PreciseStep::distance | ( | const CartesianState & | a, | |
const CartesianState & | b | |||
) | const |
Definition at line 63 of file RK4PreciseStep.cc.
References muonGeometry::mag(), CartesianState::momentum(), and CartesianState::position().
Referenced by stepWithAccuracy().
00064 { 00065 return (a.position() - b.position()).mag() + (a.momentum() - b.momentum()).mag() / b.momentum().mag(); 00066 }
CartesianState RK4PreciseStep::operator() | ( | const CartesianState & | start, | |
const RKCartesianDerivative & | deriv, | |||
double | step, | |||
double | eps | |||
) | const |
Definition at line 7 of file RK4PreciseStep.cc.
References GenMuonPlsPt100GeV_cfg::cout, lat::endl(), max, min, CartesianState::position(), funct::pow(), stepWithAccuracy(), and verbose().
00009 { 00010 const double Safety = 0.9; 00011 double remainigStep = step; 00012 double stepSize = step; 00013 CartesianState currentStart = start; 00014 int nsteps = 0; 00015 std::pair<CartesianState, double> tryStep; 00016 00017 do { 00018 tryStep = stepWithAccuracy( currentStart, deriv, stepSize); 00019 nsteps++; 00020 if (tryStep.second <eps) { 00021 if (remainigStep - stepSize < eps/2) { 00022 if (verbose()) std::cout << "Accuracy reached, and full step taken in " 00023 << nsteps << " steps" << std::endl; 00024 return tryStep.first; // we are there 00025 } 00026 else { 00027 remainigStep -= stepSize; 00028 // increase step size 00029 double factor = std::min( Safety * pow( fabs(eps/tryStep.second),0.2), 4.); 00030 stepSize = std::min( stepSize*factor, remainigStep); 00031 currentStart = tryStep.first; 00032 if (verbose()) std::cout << "Accuracy reached, but " << remainigStep 00033 << " remain after " << nsteps << " steps. Step size increased by " 00034 << factor << " to " << stepSize << std::endl; 00035 } 00036 } 00037 else { 00038 // decrease step size 00039 double factor = std::max( Safety * pow( fabs(eps/tryStep.second),0.25), 0.1); 00040 stepSize *= factor; 00041 if (verbose()) std::cout << "Accuracy not yet reached: delta = " << tryStep.second 00042 << ", step reduced by " << factor << " to " << stepSize 00043 << ", (R,z)= " << currentStart.position().perp() 00044 << ", " << currentStart.position().z() << std::endl; 00045 } 00046 } while (remainigStep > eps/2); 00047 00048 return tryStep.first; 00049 }
std::pair< CartesianState, double > RK4PreciseStep::stepWithAccuracy | ( | const CartesianState & | start, | |
const RKCartesianDerivative & | deriv, | |||
double | step | |||
) | const |
Definition at line 52 of file RK4PreciseStep.cc.
References diff, and distance().
Referenced by operator()().
00054 { 00055 RK4OneStep solver; 00056 CartesianState one(solver(start, deriv, step)); 00057 CartesianState firstHalf(solver(start, deriv, step/2)); 00058 CartesianState secondHalf(solver(firstHalf, deriv, step/2)); 00059 double diff = distance(one, secondHalf); 00060 return std::pair<CartesianState, double>(secondHalf,diff); 00061 }
Definition at line 68 of file RK4PreciseStep.cc.
Referenced by operator()().
00069 { 00070 // static bool verb = SimpleConfigurable<bool>(false,"RK4PreciseStep:verbose").value(); 00071 00072 static bool verb = true; 00073 return verb; 00074 }