TrackKinematicStatePropagator Class Reference

#include <TrackKinematicStatePropagator.h>

Inheritance diagram for TrackKinematicStatePropagator:

Public Member Functions
virtual KinematicStatePropagator *	clone () const
virtual KinematicState	propagateToTheTransversePCA (const KinematicState &state, const GlobalPoint &referencePoint) const
	TrackKinematicStatePropagator ()
virtual bool	willPropagateToTheTransversePCA (const KinematicState &state, const GlobalPoint &point) const
virtual	~TrackKinematicStatePropagator ()
Public Member Functions inherited from KinematicStatePropagator
	KinematicStatePropagator ()
virtual	~KinematicStatePropagator ()

Private Types
typedef Point3DBase< double, GlobalTag >	GlobalPointDouble
typedef Vector3DBase< double, GlobalTag >	GlobalVectorDouble

Private Member Functions
virtual KinematicState	propagateToTheTransversePCACharged (const KinematicState &state, const GlobalPoint &referencePoint) const
virtual KinematicState	propagateToTheTransversePCANeutral (const KinematicState &state, const GlobalPoint &referencePoint) const

Detailed Description

Propagator for TransientTrack based KinematicStates. Does not include the material.

Definition at line 18 of file TrackKinematicStatePropagator.h.

Member Typedef Documentation

typedef Point3DBase< double, GlobalTag> TrackKinematicStatePropagator::GlobalPointDouble

private

Definition at line 53 of file TrackKinematicStatePropagator.h.

typedef Vector3DBase< double, GlobalTag> TrackKinematicStatePropagator::GlobalVectorDouble

private

Definition at line 54 of file TrackKinematicStatePropagator.h.

Constructor & Destructor Documentation

TrackKinematicStatePropagator::TrackKinematicStatePropagator ( )

inline

Definition at line 22 of file TrackKinematicStatePropagator.h.

Referenced by clone().

{}

virtual TrackKinematicStatePropagator::~TrackKinematicStatePropagator ( )

inlinevirtual

Definition at line 24 of file TrackKinematicStatePropagator.h.

References point, propagateToTheTransversePCA(), and willPropagateToTheTransversePCA().

{}

Member Function Documentation

virtual KinematicStatePropagator* TrackKinematicStatePropagator::clone ( void  ) const

inlinevirtual

Clone method reimplemented from abstract class

Implements KinematicStatePropagator.

Definition at line 40 of file TrackKinematicStatePropagator.h.

References propagateToTheTransversePCACharged(), propagateToTheTransversePCANeutral(), and TrackKinematicStatePropagator().

  {return new TrackKinematicStatePropagator(*this);}

KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCA ( const KinematicState &  state, const GlobalPoint &  referencePoint  ) const

virtual

Propagation to the point of closest approach in transverse plane to the given point

Implements KinematicStatePropagator.

Definition at line 13 of file TrackKinematicStatePropagator.cc.

References anyDirection, Plane::build(), FreeTrajectoryState::charge(), kappa, GlobalTrajectoryParameters::magneticFieldInInverseGeV(), FreeTrajectoryState::momentum(), FreeTrajectoryState::parameters(), KinematicState::particleCharge(), point, FreeTrajectoryState::position(), makeMuonMisalignmentScenario::rot, FreeTrajectoryState::transverseCurvature(), Vector3DBase< T, FrameTag >::unit(), X, PV3DBase< T, PVType, FrameType >::x(), DOFs::Y, PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by ~TrackKinematicStatePropagator().

{
 if( state.particleCharge() == 0. ) {
    return propagateToTheTransversePCANeutral(state, referencePoint);
  } else {
    return propagateToTheTransversePCACharged(state, referencePoint);
  }
}

KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCACharged ( const KinematicState &  state, const GlobalPoint &  referencePoint  ) const

privatevirtual

Internal private methods, distinguishing between the propagation of neutrals and propagation of cahrged tracks.

Definition at line 74 of file TrackKinematicStatePropagator.cc.

References FreeTrajectoryState::curvilinearError(), HelixBarrelPlaneCrossingByCircle::direction(), KinematicState::freeTrajectoryState(), hcaldqm::quantity::fState, KinematicState::globalMomentum(), mps_fire::i, JacobianCartesianToCurvilinear::jacobian(), JacobianCurvilinearToCartesian::jacobian(), AnalyticalCurvilinearJacobian::jacobian(), KinematicState::kinematicParametersError(), LogDebug, PV3DBase< T, PVType, FrameType >::mag(), KinematicState::magneticField(), ResonanceBuilder::mass, KinematicState::mass(), KinematicParametersError::matrix(), CurvilinearTrajectoryError::matrix(), KinematicState::particleCharge(), HelixBarrelPlaneCrossingByCircle::pathLength(), HelixBarrelPlaneCrossingByCircle::position(), propagateToTheTransversePCANeutral(), alignCSCRings::s, mathSSE::sqrt(), KinematicState::trajectoryParameters(), PV3DBase< T, PVType, FrameType >::x(), Basic3DVector< T >::x(), Basic3DVector< T >::y(), and Basic3DVector< T >::z().

Referenced by clone(), and willPropagateToTheTransversePCA().

                                                                       {
  //first use the existing FTS propagator to obtain parameters at PCA
  //in transverse plane to the given point
  
  //final parameters and covariance
  
  //initial parameters as class and vectors:  
  //making a free trajectory state and looking 
  //for helix barrel plane crossing  
  FreeTrajectoryState const & fState = state.freeTrajectoryState();     
  GlobalPoint iP = referencePoint;                    
  std::pair<HelixBarrelPlaneCrossingByCircle, BoundPlane::BoundPlanePointer> cros = planeCrossing(fState,iP);      
  
  HelixBarrelPlaneCrossingByCircle planeCrossing = cros.first; 
  BoundPlane::BoundPlanePointer plane = cros.second;
  std::pair<bool,double> propResult = planeCrossing.pathLength(*plane);
  if ( !propResult.first ) {
    LogDebug("RecoVertex/TrackKinematicStatePropagator") 
    << "Propagation failed! State is invalid\n";
    return  KinematicState();
  }
  double s = propResult.second;
    
  
  GlobalTrajectoryParameters const & inPar = state.trajectoryParameters(); 
  ParticleMass mass = state.mass();                           
  GlobalVector inMom = state.globalMomentum();                            
  
  HelixPlaneCrossing::PositionType xGen = planeCrossing.position(s);
  GlobalPoint nPosition(xGen.x(),xGen.y(),xGen.z());
  HelixPlaneCrossing::DirectionType pGen = planeCrossing.direction(s);
  pGen *= inMom.mag()/pGen.mag();
  GlobalVector nMomentum(pGen.x(),pGen.y(),pGen.z()); 
  AlgebraicVector7 par;
  AlgebraicSymMatrix77 cov;
  par(0) = nPosition.x();
  par(1) = nPosition.y();
  par(2) = nPosition.z();
  par(3) = nMomentum.x();
  par(4) = nMomentum.y();
  par(5) = nMomentum.z();
  par(6) = mass;
  
  //covariance matrix business  
  //elements of 7x7 covariance matrix responcible for the mass and
  //mass - momentum projections corellations do change under such a transformation:
  //special Jacobian needed  
  GlobalTrajectoryParameters fPar(nPosition, nMomentum, state.particleCharge(),
                state.magneticField());
                                                      
  // check if correlation are present between mass and others
  bool thereIsNoCorr=true;

  for (auto i=0; i<6; ++i) thereIsNoCorr &= (0==state.kinematicParametersError().matrix()(i,6));

  if (thereIsNoCorr) {
    //  easy life
    AnalyticalCurvilinearJacobian prop(inPar,nPosition,nMomentum,s);
    AlgebraicSymMatrix55 cov2 = ROOT::Math::Similarity(prop.jacobian(), fState.curvilinearError().matrix());
    FreeTrajectoryState fts(fPar,  CurvilinearTrajectoryError(cov2));

    return  KinematicState(fts, state.mass(), std::sqrt(state.kinematicParametersError().matrix()(6,6)));

    //KinematicState kRes(fts, state.mass(), std::sqrt(state.kinematicParametersError().matrix()(6,6)));
    //std::cout << "\n\ncart from final Kstate\n" << kRes.kinematicParametersError().matrix() << std::endl;
    // std::cout << "curv from final K\n" << kRes.freeTrajectoryState().curvilinearError().matrix() << std::endl;

  } else {
  

    JacobianCartesianToCurvilinear cart2curv(inPar);
    JacobianCurvilinearToCartesian curv2cart(fPar);
    
    
    AlgebraicMatrix67 ca2cu;
    AlgebraicMatrix76 cu2ca;
    ca2cu.Place_at(cart2curv.jacobian(),0,0);
    cu2ca.Place_at(curv2cart.jacobian(),0,0);
    ca2cu(5,6) = 1;  
    cu2ca(6,5) = 1;
    
    //now both transformation jacobians: cartesian to curvilinear and back are done
    //We transform matrix to curv frame, then propagate it and translate it back to
    //cartesian frame.  
    AlgebraicSymMatrix66 cov1 = ROOT::Math::Similarity(ca2cu, state.kinematicParametersError().matrix());
    
    /*
      std::cout << "\n\ncurv from Kstate\n" << cov1 << std::endl;
      std::cout << "curv from fts\n" << fState.curvilinearError().matrix() << std::endl;
    */
    
    
    //propagation jacobian
    AnalyticalCurvilinearJacobian prop(inPar,nPosition,nMomentum,s);
    AlgebraicMatrix66 pr;
    pr(5,5) = 1;
    pr.Place_at(prop.jacobian(),0,0);

    //transportation
    AlgebraicSymMatrix66 cov2 = ROOT::Math::Similarity(pr, cov1);
    
  //now geting back to 7-parametrization from curvilinear
    cov = ROOT::Math::Similarity(cu2ca, cov2);
    
    /*
      std::cout << "curv prop \n" << cov2 << std::endl;
   std::cout << "cart prop\n" << cov << std::endl;
    */
    
    //return parameters as a kiematic state  
    KinematicParameters resPar(par);
  KinematicParametersError resEr(cov);
  
  return KinematicState(resPar,resEr,state.particleCharge(), state.magneticField());
  
  /*
    KinematicState resK(resPar,resEr,state.particleCharge(), state.magneticField());
    std::cout << "\ncurv from K prop\n" << resK.freeTrajectoryState().curvilinearError().matrix() << std::endl;
    return resK;
  */
  }

 }

KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCANeutral ( const KinematicState &  state, const GlobalPoint &  referencePoint  ) const

privatevirtual

Definition at line 199 of file TrackKinematicStatePropagator.cc.

References funct::cos(), KinematicState::freeTrajectoryState(), hcaldqm::quantity::fState, JacobianCurvilinearToCartesian::jacobian(), KinematicState::kinematicParametersError(), KinematicState::magneticField(), KinematicState::mass(), KinematicParametersError::matrix(), FreeTrajectoryState::momentum(), KinematicState::particleCharge(), PV3DBase< T, PVType, FrameType >::phi(), FreeTrajectoryState::position(), alignCSCRings::s, funct::sin(), funct::tan(), PV3DBase< T, PVType, FrameType >::theta(), theta(), KinematicState::trajectoryParameters(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by clone(), and propagateToTheTransversePCACharged().

{
//new parameters vector and covariance:
 AlgebraicVector7 par;
 AlgebraicSymMatrix77 cov;
 
 //AlgebraicVector7 inStatePar = state.kinematicParameters().vector();
 GlobalTrajectoryParameters const & inPar = state.trajectoryParameters(); 
 
 //first making a free trajectory state and propagating it according
 //to the algorithm provided by Thomas Speer and Wolfgang Adam 
 FreeTrajectoryState const &  fState = state.freeTrajectoryState();
  
 GlobalPoint xvecOrig = fState.position();
 double phi = fState.momentum().phi();
 double theta = fState.momentum().theta();
 double xOrig = xvecOrig.x();
 double yOrig = xvecOrig.y();
 double zOrig = xvecOrig.z();
 double xR = referencePoint.x();
 double yR = referencePoint.y();

 double s2D = (xR - xOrig)  * cos(phi) + (yR - yOrig)  * sin(phi);
 double s = s2D / sin(theta);
 double xGen = xOrig + s2D*cos(phi);
 double yGen = yOrig + s2D*sin(phi);
 double zGen = zOrig + s2D/tan(theta);
 GlobalPoint xPerigee = GlobalPoint(xGen, yGen, zGen);

//new parameters
 GlobalVector pPerigee = fState.momentum();
 par(0) = xPerigee.x();
 par(1) = xPerigee.y(); 
 par(2) = xPerigee.z(); 
 par(3) = pPerigee.x(); 
 par(4) = pPerigee.y(); 
 par(5) = pPerigee.z(); 
 // par(6) = inStatePar(7);
 par(6) = state.mass();

//covariance matrix business:
//everything lake it was before: jacobains are smart enouhg to 
//distinguish between neutral and charged states themselves

 GlobalTrajectoryParameters fPar(xPerigee, pPerigee, state.particleCharge(),
                state.magneticField());

 JacobianCartesianToCurvilinear cart2curv(inPar);
 JacobianCurvilinearToCartesian curv2cart(fPar);
  
  AlgebraicMatrix67 ca2cu;
  AlgebraicMatrix76 cu2ca;
  ca2cu.Place_at(cart2curv.jacobian(),0,0);
  cu2ca.Place_at(curv2cart.jacobian(),0,0);
  ca2cu(5,6) = 1;  
  cu2ca(6,5) = 1;
  
  //now both transformation jacobians: cartesian to curvilinear and back are done
  //We transform matrix to curv frame, then propagate it and translate it back to
  //cartesian frame.  
  AlgebraicSymMatrix66 cov1 = ROOT::Math::Similarity(ca2cu, state.kinematicParametersError().matrix());
  
  //propagation jacobian
  AnalyticalCurvilinearJacobian prop(inPar,xPerigee,pPerigee,s);
  AlgebraicMatrix66 pr;
  pr(5,5) = 1;
  pr.Place_at(prop.jacobian(),0,0);
  
  //transportation
  AlgebraicSymMatrix66 cov2 = ROOT::Math::Similarity(pr, cov1);
  
  //now geting back to 7-parametrization from curvilinear
  cov = ROOT::Math::Similarity(cu2ca, cov2);
  
  // FreeTrajectoryState fts(fPar);
  
  //return parameters as a kiematic state  
  KinematicParameters resPar(par);
  KinematicParametersError resEr(cov);
  return  KinematicState(resPar,resEr,state.particleCharge(), state.magneticField()); 

  //return  KinematicState(fts,state.mass(), cov(6,6));

}

bool TrackKinematicStatePropagator::willPropagateToTheTransversePCA ( const KinematicState &  state, const GlobalPoint &  point  ) const

virtual

Reimplemented from KinematicStatePropagator.

Definition at line 58 of file TrackKinematicStatePropagator.cc.

References KinematicState::freeTrajectoryState(), hcaldqm::quantity::fState, KinematicState::particleCharge(), HelixBarrelPlaneCrossingByCircle::pathLength(), and propagateToTheTransversePCACharged().

Referenced by KinematicParticleVertexFitter::fit(), and ~TrackKinematicStatePropagator().

                                                                                                                               {
  if( state.particleCharge() == 0. ) return true;

  // copied from below...
  FreeTrajectoryState const & fState = state.freeTrajectoryState();     
  std::pair<HelixBarrelPlaneCrossingByCircle, BoundPlane::BoundPlanePointer> cros = planeCrossing(fState,point);       
  
  HelixBarrelPlaneCrossingByCircle planeCrossing = cros.first; 
  BoundPlane::BoundPlanePointer plane = cros.second;
  std::pair<bool,double> propResult = planeCrossing.pathLength(*plane);
  return propResult.first;

}