|
|
|
#include <RKPropagatorInS.h>
Definition at line 20 of file RKPropagatorInS.h.
typedef std::pair<TrajectoryStateOnSurface,double> RKPropagatorInS::TsosWP [private] |
Definition at line 61 of file RKPropagatorInS.h.
| RKPropagatorInS::RKPropagatorInS | ( | const MagVolume & | vol, |
| PropagationDirection | dir = alongMomentum, |
||
| double | tolerance = 5.e-5 |
||
| ) | [inline, explicit] |
Definition at line 27 of file RKPropagatorInS.h.
Referenced by clone().
:
Propagator(dir), theVolume( &vol), theTolerance( tolerance) {}
| RKPropagatorInS::~RKPropagatorInS | ( | ) | [inline] |
Definition at line 31 of file RKPropagatorInS.h.
{}
| Propagator * RKPropagatorInS::clone | ( | void | ) | const [virtual] |
Implements Propagator.
Definition at line 345 of file RKPropagatorInS.cc.
References RKPropagatorInS().
{
return new RKPropagatorInS(*this);
}
| RKLocalFieldProvider RKPropagatorInS::fieldProvider | ( | ) | const [private] |
Definition at line 358 of file RKPropagatorInS.cc.
References theVolume.
Referenced by propagateParametersOnCylinder(), and propagateParametersOnPlane().
{
return RKLocalFieldProvider( *theVolume);
}
| RKLocalFieldProvider RKPropagatorInS::fieldProvider | ( | const Cylinder & | cyl | ) | const [private] |
Definition at line 363 of file RKPropagatorInS.cc.
References theVolume.
{
return RKLocalFieldProvider( *theVolume, cyl);
}
| GlobalVector RKPropagatorInS::globalMomentum | ( | const Basic3DVector< double > & | mom | ) | const [private] |
Definition at line 392 of file RKPropagatorInS.cc.
References theVolume, and GloballyPositioned< T >::toGlobal().
Referenced by gtpFromVolumeLocal().
{
if (theVolume != 0) return theVolume->toGlobal( LocalVector(mom));
else return GlobalVector(mom);
}
| GlobalPoint RKPropagatorInS::globalPosition | ( | const Basic3DVector< double > & | pos | ) | const [private] |
Definition at line 386 of file RKPropagatorInS.cc.
References theVolume, and GloballyPositioned< T >::toGlobal().
Referenced by gtpFromVolumeLocal(), and propagateParametersOnPlane().
{
if (theVolume != 0) return theVolume->toGlobal( LocalPoint(pos));
else return GlobalPoint(pos);
}
| GlobalTrajectoryParameters RKPropagatorInS::gtpFromLocal | ( | const Basic3DVector< double > & | lpos, |
| const Basic3DVector< double > & | lmom, | ||
| TrackCharge | ch, | ||
| const Surface & | surf | ||
| ) | const [private] |
Definition at line 350 of file RKPropagatorInS.cc.
References theVolume, and Surface::toGlobal().
Referenced by propagateParametersOnCylinder().
{
return GlobalTrajectoryParameters( surf.toGlobal( LocalPoint( lpos)),
surf.toGlobal( LocalVector( lmom)), ch, theVolume);
}
| GlobalTrajectoryParameters RKPropagatorInS::gtpFromVolumeLocal | ( | const CartesianStateAdaptor & | state, |
| TrackCharge | charge | ||
| ) | const [private] |
Definition at line 400 of file RKPropagatorInS.cc.
References globalMomentum(), globalPosition(), CartesianStateAdaptor::momentum(), CartesianStateAdaptor::position(), and theVolume.
Referenced by propagateParametersOnPlane().
{
return GlobalTrajectoryParameters( globalPosition(state.position()),
globalMomentum(state.momentum()),
charge, theVolume);
}
| PropagationDirection RKPropagatorInS::invertDirection | ( | PropagationDirection | dir | ) | const [private] |
Definition at line 368 of file RKPropagatorInS.cc.
References alongMomentum, anyDirection, dir, and oppositeToMomentum.
Referenced by propagateParametersOnCylinder(), and propagateParametersOnPlane().
{
if (dir == anyDirection) return dir;
return ( dir == alongMomentum ? oppositeToMomentum : alongMomentum);
}
| virtual const MagneticField* RKPropagatorInS::magneticField | ( | ) | const [inline, virtual] |
Implements Propagator.
Definition at line 57 of file RKPropagatorInS.h.
References theVolume.
Referenced by NavPropagator::propagateInVolume().
{return theVolume;}
| TrajectoryStateOnSurface RKPropagatorInS::propagate | ( | const FreeTrajectoryState & | fts, |
| const Cylinder & | cyl | ||
| ) | const [virtual] |
Implements Propagator.
Definition at line 340 of file RKPropagatorInS.cc.
References propagateWithPath().
{
return propagateWithPath( fts, cyl).first;
}
| TrajectoryStateOnSurface RKPropagatorInS::propagate | ( | const TrajectoryStateOnSurface & | ts, |
| const Plane & | plane | ||
| ) | const [inline, virtual] |
Reimplemented from Propagator.
Definition at line 48 of file RKPropagatorInS.h.
References first, TrajectoryStateOnSurface::freeState(), and propagateWithPath().
{
return propagateWithPath( *ts.freeState(),plane).first;
}
| TrajectoryStateOnSurface RKPropagatorInS::propagate | ( | const FreeTrajectoryState & | fts, |
| const Plane & | plane | ||
| ) | const [virtual] |
Implements Propagator.
Definition at line 334 of file RKPropagatorInS.cc.
References propagateWithPath().
{
return propagateWithPath( fts, plane).first;
}
| GlobalParametersWithPath RKPropagatorInS::propagateParametersOnCylinder | ( | const FreeTrajectoryState & | ts, |
| const Cylinder & | cyl | ||
| ) | const [private] |
Definition at line 200 of file RKPropagatorInS.cc.
References abs, FreeTrajectoryState::charge(), fieldProvider(), gtpFromLocal(), invertDirection(), likely, LogDebug, FreeTrajectoryState::momentum(), CartesianStateAdaptor::momentum(), path(), StraightLineBarrelCylinderCrossing::pathLength(), StraightLineCylinderCrossing::pathLength(), perp(), pos, CartesianStateAdaptor::position(), StraightLineBarrelCylinderCrossing::position(), GloballyPositioned< T >::position(), FreeTrajectoryState::position(), Propagator::propagationDirection(), Cylinder::radius(), rho, CartesianStateAdaptor::rkstate(), asciidump::s, theTolerance, theVolume, GloballyPositioned< T >::toLocal(), FreeTrajectoryState::transverseCurvature(), unlikely, PV3DBase< T, PVType, FrameType >::x(), x, and PV3DBase< T, PVType, FrameType >::y().
Referenced by propagateWithPath().
{
typedef RKAdaptiveSolver<double,RKOneCashKarpStep, 6> Solver;
typedef Solver::Vector RKVector;
GlobalPoint sp = cyl.position();
if unlikely(sp.x()!=0. || sp.y()!=0.) {
throw PropagationException("Cannot propagate to an arbitrary cylinder");
}
GlobalPoint gpos( ts.position());
GlobalVector gmom( ts.momentum());
LocalPoint pos(cyl.toLocal(gpos));
LocalVector mom(cyl.toLocal(gmom));
double startR = cyl.radius() - pos.perp();
// LogDebug("RKPropagatorInS") << "RKPropagatorInS: starting from FTS " << ts ;
// (transverse) curvature
double rho = ts.transverseCurvature();
//
// Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
// difference in transversal position at 10m.
//
if unlikely( fabs(rho)<1.e-10 ) {
//
// Instantiate auxiliary object for finding intersection.
// Frame-independant point and vector are created explicitely to
// avoid confusing gcc (refuses to compile with temporary objects
// in the constructor).
//
StraightLineBarrelCylinderCrossing cylCrossing(gpos,gmom,propagationDirection());
//
// get solution
//
std::pair<bool,double> propResult = cylCrossing.pathLength(cyl);
if likely( propResult.first && theVolume !=0) {
double s = propResult.second;
// point (reconverted to GlobalPoint)
GlobalPoint x (cylCrossing.position(s));
GlobalTrajectoryParameters res( x, gmom, ts.charge(), theVolume);
LogDebug("RKPropagatorInS") << "Straight line propagation to cylinder succeeded !!";
return GlobalParametersWithPath( res, s);
}
//do someting
edm::LogError("RKPropagatorInS") <<"Straight line propagation to cylinder failed !!";
return GlobalParametersWithPath( );
}
RKLocalFieldProvider field( fieldProvider(cyl));
// StraightLineCylinderCrossing pathLength( pos, mom, propagationDirection());
CartesianLorentzForce deriv(field, ts.charge());
RKCartesianDistance dist;
double eps = theTolerance;
Solver solver;
double stot = 0;
PropagationDirection currentDirection = propagationDirection();
RKVector start( CartesianStateAdaptor::rkstate( pos.basicVector(), mom.basicVector()));
int safeGuard = 0;
while (safeGuard++<100) {
CartesianStateAdaptor startState(start);
StraightLineCylinderCrossing pathLength( LocalPoint(startState.position()),
LocalVector(startState.momentum()),
currentDirection, eps);
std::pair<bool,double> path = pathLength.pathLength( cyl);
if unlikely(!path.first) {
LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path length calculation to cylinder failed!"
<< "Radius " << cyl.radius() << " pos.perp() " << LocalPoint(startState.position()).perp() ;
return GlobalParametersWithPath();
}
LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path lenght to cylinder is " << path.second
<< " from point (R,z) " << startState.position().perp()
<< ", " << startState.position().z()
<< " to R " << cyl.radius()
;
double sstep = path.second;
if unlikely( std::abs(sstep) < eps) {
LogDebug("RKPropagatorInS") << "accuracy not reached, but pathLength calculation says we are there! "
<< path.second ;
GlobalTrajectoryParameters res( gtpFromLocal( startState.position(),
startState.momentum(),
ts.charge(), cyl));
return GlobalParametersWithPath( res, stot);
}
LogDebug("RKPropagatorInS") << "RKPropagatorInS: Solving for " << sstep
<< " current distance to cylinder is " << startR ;
RKVector rkresult = solver( 0, start, sstep, deriv, dist, eps);
stot += sstep;
CartesianStateAdaptor cur( rkresult);
double remainingR = cyl.radius() - cur.position().perp();
if ( fabs(remainingR) < eps) {
LogDebug("RKPropagatorInS") << "Accuracy reached! " << remainingR ;
GlobalTrajectoryParameters res( gtpFromLocal( cur.position(),
cur.momentum(),
ts.charge(), cyl));
return GlobalParametersWithPath( res, stot);
}
start = rkresult;
if (remainingR * startR > 0) {
LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in same direction again "
<< remainingR ;
}
else {
LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in opposite direction "
<< remainingR ;
currentDirection = invertDirection( currentDirection);
}
startR = remainingR;
}
edm::LogError("FailedPropagation") << " too many iterations trying to reach cylinder ";
return GlobalParametersWithPath();
}
| GlobalParametersWithPath RKPropagatorInS::propagateParametersOnPlane | ( | const FreeTrajectoryState & | ts, |
| const Plane & | plane | ||
| ) | const [private] |
Definition at line 53 of file RKPropagatorInS.cc.
References abs, PV3DBase< T, PVType, FrameType >::basicVector(), FreeTrajectoryState::charge(), cross(), dir, fieldProvider(), RKLocalFieldProvider::frame(), globalPosition(), gtpFromVolumeLocal(), RKLocalFieldProvider::inTesla(), invertDirection(), likely, Plane::localZ(), LogDebug, FreeTrajectoryState::momentum(), CartesianStateAdaptor::momentum(), path(), StraightLinePlaneCrossing::pathLength(), pos, CartesianStateAdaptor::position(), StraightLinePlaneCrossing::position(), GloballyPositioned< T >::position(), FreeTrajectoryState::position(), Propagator::propagationDirection(), rho, rkMomentum(), rkPosition(), CartesianStateAdaptor::rkstate(), asciidump::s, theTolerance, theVolume, GloballyPositioned< T >::toGlobal(), GloballyPositioned< T >::toLocal(), FrameChanger::transformPlane(), FreeTrajectoryState::transverseCurvature(), unlikely, and x.
Referenced by propagateWithPath().
{
GlobalPoint gpos( ts.position());
GlobalVector gmom( ts.momentum());
double startZ = plane.localZ(gpos);
// (transverse) curvature
double rho = ts.transverseCurvature();
//
// Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
// difference in transversal position at 10m.
//
if unlikely( fabs(rho)<1.e-10 ) {
//
// Instantiate auxiliary object for finding intersection.
// Frame-independant point and vector are created explicitely to
// avoid confusing gcc (refuses to compile with temporary objects
// in the constructor).
//
LogDebug("RKPropagatorInS")<<" startZ = "<<startZ;
if unlikely(fabs(startZ) < 1e-5){
LogDebug("RKPropagatorInS")<< "Propagation is not performed: state is already on final surface.";
GlobalTrajectoryParameters res( gpos, gmom, ts.charge(),
theVolume);
return GlobalParametersWithPath( res, 0.0);
}
StraightLinePlaneCrossing::PositionType pos(gpos);
StraightLinePlaneCrossing::DirectionType dir(gmom);
StraightLinePlaneCrossing planeCrossing(pos,dir, propagationDirection());
//
// get solution
//
std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
if likely( propResult.first && theVolume !=0) {
double s = propResult.second;
// point (reconverted to GlobalPoint)
GlobalPoint x (planeCrossing.position(s));
GlobalTrajectoryParameters res( x, gmom, ts.charge(), theVolume);
return GlobalParametersWithPath( res, s);
}
//do someting
LogDebug("RKPropagatorInS")<< "Straight line propgation to plane failed !!";
return GlobalParametersWithPath( );
}
#ifdef EDM_LM_DEBUG
if (theVolume != 0) {
LogDebug("RKPropagatorInS") << "RKPropagatorInS: starting prop to plane in volume with pos " << theVolume->position()
<< " Z axis " << theVolume->toGlobal( LocalVector(0,0,1)) ;
LogDebug("RKPropagatorInS") << "The starting position is " << ts.position() << " (global) "
<< theVolume->toLocal(ts.position()) << " (local) " ;
FrameChanger changer;
FrameChanger::PlanePtr localPlane = changer.transformPlane( plane, *theVolume);
LogDebug("RKPropagatorInS") << "The plane position is " << plane.position() << " (global) "
<< localPlane->position() << " (local) " ;
LogDebug("RKPropagatorInS") << "The initial distance to plane is " << plane.localZ( ts.position()) ;
StraightLinePlaneCrossing cross( ts.position().basicVector(), ts.momentum().basicVector());
std::pair<bool,double> res3 = cross.pathLength(plane);
LogDebug("RKPropagatorInS") << "straight line distance " << res3.first << " " << res3.second ;
}
#endif
typedef RKAdaptiveSolver<double,RKOneCashKarpStep, 6> Solver;
typedef Solver::Vector RKVector;
RKLocalFieldProvider field( fieldProvider());
PathToPlane2Order pathLength( field, &field.frame());
CartesianLorentzForce deriv(field, ts.charge());
RKCartesianDistance dist;
double eps = theTolerance;
Solver solver;
double stot = 0;
PropagationDirection currentDirection = propagationDirection();
// in magVolume frame
RKVector start( CartesianStateAdaptor::rkstate( rkPosition(gpos), rkMomentum(gmom)));
int safeGuard = 0;
while (safeGuard++<100) {
CartesianStateAdaptor startState(start);
std::pair<bool,double> path = pathLength( plane, startState.position(),
startState.momentum(),
(double) ts.charge(), currentDirection);
if unlikely(!path.first) {
LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path length calculation to plane failed!"
<< "...distance to plane " << plane.localZ( globalPosition(startState.position()))
<< "...Local starting position in volume " << startState.position()
<< "...Magnetic field " << field.inTesla( startState.position()) ;
return GlobalParametersWithPath();
}
LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path lenght to plane is " << path.second ;
double sstep = path.second;
if unlikely( std::abs(sstep) < eps) {
LogDebug("RKPropagatorInS") << "On-surface accuracy not reached, but pathLength calculation says we are there! "
<< "path " << path.second << " distance to plane is " << startZ ;
GlobalTrajectoryParameters res( gtpFromVolumeLocal( startState, ts.charge()));
return GlobalParametersWithPath( res, stot);
}
LogDebug("RKPropagatorInS") << "RKPropagatorInS: Solving for " << sstep
<< " current distance to plane is " << startZ ;
RKVector rkresult = solver( 0, start, sstep, deriv, dist, eps);
stot += sstep;
CartesianStateAdaptor cur( rkresult);
double remainingZ = plane.localZ( globalPosition(cur.position()));
if ( fabs(remainingZ) < eps) {
LogDebug("RKPropagatorInS") << "On-surface accuracy reached! " << remainingZ ;
GlobalTrajectoryParameters res( gtpFromVolumeLocal( cur, ts.charge()));
return GlobalParametersWithPath( res, stot);
}
start = rkresult;
if (remainingZ * startZ > 0) {
LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in same direction again "
<< remainingZ ;
}
else {
LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in opposite direction "
<< remainingZ ;
currentDirection = invertDirection( currentDirection);
}
startZ = remainingZ;
}
edm::LogError("FailedPropagation") << " too many iterations trying to reach plane ";
return GlobalParametersWithPath();
}
| std::pair< TrajectoryStateOnSurface, double > RKPropagatorInS::propagateWithPath | ( | const FreeTrajectoryState & | fts, |
| const Cylinder & | cyl | ||
| ) | const [virtual] |
Implements Propagator.
Definition at line 40 of file RKPropagatorInS.cc.
References SurfaceSideDefinition::afterSurface, alongMomentum, SurfaceSideDefinition::beforeSurface, GlobalParametersWithPath::parameters(), propagateParametersOnCylinder(), Propagator::propagationDirection(), GlobalParametersWithPath::s(), and unlikely.
{
GlobalParametersWithPath gp = propagateParametersOnCylinder( fts, cyl);
if unlikely(!gp) return TsosWP(TrajectoryStateOnSurface(),0.);
SurfaceSideDefinition::SurfaceSide side = PropagationDirectionFromPath()(gp.s(),propagationDirection())==alongMomentum
? SurfaceSideDefinition::beforeSurface : SurfaceSideDefinition::afterSurface;
AnalyticalErrorPropagation errorprop;
return errorprop( fts, cyl, side, gp.parameters(),gp.s());
}
| std::pair< TrajectoryStateOnSurface, double > RKPropagatorInS::propagateWithPath | ( | const FreeTrajectoryState & | fts, |
| const Plane & | plane | ||
| ) | const [virtual] |
Implements Propagator.
Definition at line 27 of file RKPropagatorInS.cc.
References SurfaceSideDefinition::afterSurface, alongMomentum, SurfaceSideDefinition::beforeSurface, GlobalParametersWithPath::parameters(), propagateParametersOnPlane(), Propagator::propagationDirection(), GlobalParametersWithPath::s(), and unlikely.
Referenced by propagate(), and NavPropagator::propagateInVolume().
{
GlobalParametersWithPath gp = propagateParametersOnPlane( fts, plane);
if unlikely(!gp) return TsosWP(TrajectoryStateOnSurface(),0.);
SurfaceSideDefinition::SurfaceSide side = PropagationDirectionFromPath()(gp.s(),propagationDirection())==alongMomentum
? SurfaceSideDefinition::beforeSurface : SurfaceSideDefinition::afterSurface;
AnalyticalErrorPropagation errorprop;
return errorprop( fts, plane, side, gp.parameters(),gp.s());
}
| Basic3DVector< double > RKPropagatorInS::rkMomentum | ( | const GlobalVector & | mom | ) | const [private] |
Definition at line 380 of file RKPropagatorInS.cc.
References PV3DBase< T, PVType, FrameType >::basicVector(), theVolume, and GloballyPositioned< T >::toLocal().
Referenced by propagateParametersOnPlane().
{
if (theVolume != 0) return theVolume->toLocal( mom).basicVector();
else return mom.basicVector();
}
| Basic3DVector< double > RKPropagatorInS::rkPosition | ( | const GlobalPoint & | pos | ) | const [private] |
Definition at line 374 of file RKPropagatorInS.cc.
References PV3DBase< T, PVType, FrameType >::basicVector(), theVolume, and GloballyPositioned< T >::toLocal().
Referenced by propagateParametersOnPlane().
{
if (theVolume != 0) return theVolume->toLocal( pos).basicVector();
else return pos.basicVector();
}
double RKPropagatorInS::theTolerance [private] |
Definition at line 64 of file RKPropagatorInS.h.
Referenced by propagateParametersOnCylinder(), and propagateParametersOnPlane().
const MagVolume* RKPropagatorInS::theVolume [private] |
Definition at line 63 of file RKPropagatorInS.h.
Referenced by fieldProvider(), globalMomentum(), globalPosition(), gtpFromLocal(), gtpFromVolumeLocal(), magneticField(), propagateParametersOnCylinder(), propagateParametersOnPlane(), rkMomentum(), and rkPosition().
1.7.3