KinematicFitPrimitives/src/TrackKinematicStatePropagator.cc

0001 #include "RecoVertex/KinematicFitPrimitives/interface/TrackKinematicStatePropagator.h"
0002 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToCurvilinear.h"
0003 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToCartesian.h"
0004 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0006 #include "DataFormats/GeometrySurface/interface/OpenBounds.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 using namespace std;
0010
0011 KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCA(const KinematicState& state,
0012                                                                           const GlobalPoint& referencePoint) const {
0013   if (state.particleCharge() == 0.) {
0014     return propagateToTheTransversePCANeutral(state, referencePoint);
0015   } else {
0016     return propagateToTheTransversePCACharged(state, referencePoint);
0017   }
0018 }
0019
0020 namespace {
0021   inline pair<HelixBarrelPlaneCrossingByCircle, BoundPlane::BoundPlanePointer> planeCrossing(
0022       const FreeTrajectoryState& state, const GlobalPoint& point) {
0023     typedef Point3DBase<double, GlobalTag> GlobalPointDouble;
0024     typedef Vector3DBase<double, GlobalTag> GlobalVectorDouble;
0025
0026     GlobalPoint inPos = state.position();
0027     GlobalVector inMom = state.momentum();
0028     double kappa = state.transverseCurvature();
0029     auto bz = state.parameters().magneticFieldInInverseGeV(point).z();
0030     if (std::abs(bz) < 1e-6) {
0031       LogDebug("RecoVertex/TrackKinematicStatePropagator") << "planeCrossing is not possible";
0032       return {HelixBarrelPlaneCrossingByCircle(inPos, inMom, kappa), BoundPlane::BoundPlanePointer()};
0033     }
0034     double fac = state.charge() / bz;
0035
0036     GlobalVectorDouble xOrig2Centre(fac * inMom.y(), -fac * inMom.x(), 0.);
0037     GlobalVectorDouble xOrigProj(inPos.x(), inPos.y(), 0.);
0038     GlobalVectorDouble xRefProj(point.x(), point.y(), 0.);
0039     GlobalVectorDouble deltax = xRefProj - xOrigProj - xOrig2Centre;
0040     GlobalVectorDouble ndeltax = deltax.unit();
0041
0042     PropagationDirection direction = anyDirection;
0043     Surface::PositionType pos(point);
0044
0045     // Need to define plane with orientation as the ImpactPointSurface
0046     GlobalVector X(ndeltax.x(), ndeltax.y(), ndeltax.z());
0047     GlobalVector Y(0., 0., 1.);
0048     Surface::RotationType rot(X, Y);
0049     Plane::PlanePointer plane = Plane::build(pos, rot);
0050     HelixBarrelPlaneCrossingByCircle planeCrossing(HelixPlaneCrossing::PositionType(inPos.x(), inPos.y(), inPos.z()),
0051                                                    HelixPlaneCrossing::DirectionType(inMom.x(), inMom.y(), inMom.z()),
0052                                                    kappa,
0053                                                    direction);
0054     return std::pair<HelixBarrelPlaneCrossingByCircle, Plane::PlanePointer>(planeCrossing, plane);
0055   }
0056 }  // namespace
0057
0058 bool TrackKinematicStatePropagator::willPropagateToTheTransversePCA(const KinematicState& state,
0059                                                                     const GlobalPoint& point) const {
0060   if (state.particleCharge() == 0.)
0061     return true;
0062
0063   // copied from below...
0064   FreeTrajectoryState const& fState = state.freeTrajectoryState();
0065   std::pair<HelixBarrelPlaneCrossingByCircle, BoundPlane::BoundPlanePointer> cros = planeCrossing(fState, point);
0066   if (cros.second == nullptr)
0067     return false;
0068
0069   HelixBarrelPlaneCrossingByCircle planeCrossing = cros.first;
0070   BoundPlane::BoundPlanePointer plane = cros.second;
0071   std::pair<bool, double> propResult = planeCrossing.pathLength(*plane);
0072   return propResult.first;
0073 }
0074
0075 KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCACharged(
0076     const KinematicState& state, const GlobalPoint& referencePoint) const {
0077   //first use the existing FTS propagator to obtain parameters at PCA
0078   //in transverse plane to the given point
0079
0080   //final parameters and covariance
0081
0082   //initial parameters as class and vectors:
0083   //making a free trajectory state and looking
0084   //for helix barrel plane crossing
0085   FreeTrajectoryState const& fState = state.freeTrajectoryState();
0086   const GlobalPoint& iP = referencePoint;
0087   std::pair<HelixBarrelPlaneCrossingByCircle, BoundPlane::BoundPlanePointer> cros = planeCrossing(fState, iP);
0088   if (cros.second == nullptr)
0089     return KinematicState();
0090
0091   HelixBarrelPlaneCrossingByCircle planeCrossing = cros.first;
0092   BoundPlane::BoundPlanePointer plane = cros.second;
0093   std::pair<bool, double> propResult = planeCrossing.pathLength(*plane);
0094   if (!propResult.first) {
0095     LogDebug("RecoVertex/TrackKinematicStatePropagator") << "Propagation failed! State is invalid\n";
0096     return KinematicState();
0097   }
0098   double s = propResult.second;
0099
0100   GlobalTrajectoryParameters const& inPar = state.trajectoryParameters();
0101   ParticleMass mass = state.mass();
0102   GlobalVector inMom = state.globalMomentum();
0103
0104   HelixPlaneCrossing::PositionType xGen = planeCrossing.position(s);
0105   GlobalPoint nPosition(xGen.x(), xGen.y(), xGen.z());
0106   HelixPlaneCrossing::DirectionType pGen = planeCrossing.direction(s);
0107   pGen *= inMom.mag() / pGen.mag();
0108   GlobalVector nMomentum(pGen.x(), pGen.y(), pGen.z());
0109   AlgebraicVector7 par;
0110   AlgebraicSymMatrix77 cov;
0111   par(0) = nPosition.x();
0112   par(1) = nPosition.y();
0113   par(2) = nPosition.z();
0114   par(3) = nMomentum.x();
0115   par(4) = nMomentum.y();
0116   par(5) = nMomentum.z();
0117   par(6) = mass;
0118
0119   //covariance matrix business
0120   //elements of 7x7 covariance matrix responcible for the mass and
0121   //mass - momentum projections corellations do change under such a transformation:
0122   //special Jacobian needed
0123   GlobalTrajectoryParameters fPar(nPosition, nMomentum, state.particleCharge(), state.magneticField());
0124
0125   // check if correlation are present between mass and others
0126   bool thereIsNoCorr = true;
0127
0128   for (auto i = 0; i < 6; ++i)
0129     thereIsNoCorr &= (0 == state.kinematicParametersError().matrix()(i, 6));
0130
0131   if (thereIsNoCorr) {
0132     //  easy life
0133     AnalyticalCurvilinearJacobian prop(inPar, nPosition, nMomentum, s);
0134     AlgebraicSymMatrix55 cov2 = ROOT::Math::Similarity(prop.jacobian(), fState.curvilinearError().matrix());
0135     FreeTrajectoryState fts(fPar, CurvilinearTrajectoryError(cov2));
0136
0137     return KinematicState(fts, state.mass(), std::sqrt(state.kinematicParametersError().matrix()(6, 6)));
0138
0139     //KinematicState kRes(fts, state.mass(), std::sqrt(state.kinematicParametersError().matrix()(6,6)));
0140     //std::cout << "\n\ncart from final Kstate\n" << kRes.kinematicParametersError().matrix() << std::endl;
0141     // std::cout << "curv from final K\n" << kRes.freeTrajectoryState().curvilinearError().matrix() << std::endl;
0142
0143   } else {
0144     JacobianCartesianToCurvilinear cart2curv(inPar);
0145     JacobianCurvilinearToCartesian curv2cart(fPar);
0146
0147     AlgebraicMatrix67 ca2cu;
0148     AlgebraicMatrix76 cu2ca;
0149     ca2cu.Place_at(cart2curv.jacobian(), 0, 0);
0150     cu2ca.Place_at(curv2cart.jacobian(), 0, 0);
0151     ca2cu(5, 6) = 1;
0152     cu2ca(6, 5) = 1;
0153
0154     //now both transformation jacobians: cartesian to curvilinear and back are done
0155     //We transform matrix to curv frame, then propagate it and translate it back to
0156     //cartesian frame.
0157     AlgebraicSymMatrix66 cov1 = ROOT::Math::Similarity(ca2cu, state.kinematicParametersError().matrix());
0158
0159     /*
0160       std::cout << "\n\ncurv from Kstate\n" << cov1 << std::endl;
0161       std::cout << "curv from fts\n" << fState.curvilinearError().matrix() << std::endl;
0162     */
0163
0164     //propagation jacobian
0165     AnalyticalCurvilinearJacobian prop(inPar, nPosition, nMomentum, s);
0166     AlgebraicMatrix66 pr;
0167     pr(5, 5) = 1;
0168     pr.Place_at(prop.jacobian(), 0, 0);
0169
0170     //transportation
0171     AlgebraicSymMatrix66 cov2 = ROOT::Math::Similarity(pr, cov1);
0172
0173     //now geting back to 7-parametrization from curvilinear
0174     cov = ROOT::Math::Similarity(cu2ca, cov2);
0175
0176     /*
0177       std::cout << "curv prop \n" << cov2 << std::endl;
0178    std::cout << "cart prop\n" << cov << std::endl;
0179     */
0180
0181     //return parameters as a kiematic state
0182     KinematicParameters resPar(par);
0183     KinematicParametersError resEr(cov);
0184
0185     return KinematicState(resPar, resEr, state.particleCharge(), state.magneticField());
0186
0187     /*
0188     KinematicState resK(resPar,resEr,state.particleCharge(), state.magneticField());
0189     std::cout << "\ncurv from K prop\n" << resK.freeTrajectoryState().curvilinearError().matrix() << std::endl;
0190     return resK;
0191   */
0192   }
0193 }
0194
0195 KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCANeutral(
0196     const KinematicState& state, const GlobalPoint& referencePoint) const {
0197   //new parameters vector and covariance:
0198   AlgebraicVector7 par;
0199   AlgebraicSymMatrix77 cov;
0200
0201   //AlgebraicVector7 inStatePar = state.kinematicParameters().vector();
0202   GlobalTrajectoryParameters const& inPar = state.trajectoryParameters();
0203
0204   //first making a free trajectory state and propagating it according
0205   //to the algorithm provided by Thomas Speer and Wolfgang Adam
0206   FreeTrajectoryState const& fState = state.freeTrajectoryState();
0207
0208   GlobalPoint xvecOrig = fState.position();
0209   double phi = fState.momentum().phi();
0210   double theta = fState.momentum().theta();
0211   double xOrig = xvecOrig.x();
0212   double yOrig = xvecOrig.y();
0213   double zOrig = xvecOrig.z();
0214   double xR = referencePoint.x();
0215   double yR = referencePoint.y();
0216
0217   double s2D = (xR - xOrig) * cos(phi) + (yR - yOrig) * sin(phi);
0218   double s = s2D / sin(theta);
0219   double xGen = xOrig + s2D * cos(phi);
0220   double yGen = yOrig + s2D * sin(phi);
0221   double zGen = zOrig + s2D / tan(theta);
0222   GlobalPoint xPerigee = GlobalPoint(xGen, yGen, zGen);
0223
0224   //new parameters
0225   GlobalVector pPerigee = fState.momentum();
0226   par(0) = xPerigee.x();
0227   par(1) = xPerigee.y();
0228   par(2) = xPerigee.z();
0229   par(3) = pPerigee.x();
0230   par(4) = pPerigee.y();
0231   par(5) = pPerigee.z();
0232   // par(6) = inStatePar(7);
0233   par(6) = state.mass();
0234
0235   //covariance matrix business:
0236   //everything lake it was before: jacobains are smart enouhg to
0237   //distinguish between neutral and charged states themselves
0238
0239   GlobalTrajectoryParameters fPar(xPerigee, pPerigee, state.particleCharge(), state.magneticField());
0240
0241   JacobianCartesianToCurvilinear cart2curv(inPar);
0242   JacobianCurvilinearToCartesian curv2cart(fPar);
0243
0244   AlgebraicMatrix67 ca2cu;
0245   AlgebraicMatrix76 cu2ca;
0246   ca2cu.Place_at(cart2curv.jacobian(), 0, 0);
0247   cu2ca.Place_at(curv2cart.jacobian(), 0, 0);
0248   ca2cu(5, 6) = 1;
0249   cu2ca(6, 5) = 1;
0250
0251   //now both transformation jacobians: cartesian to curvilinear and back are done
0252   //We transform matrix to curv frame, then propagate it and translate it back to
0253   //cartesian frame.
0254   AlgebraicSymMatrix66 cov1 = ROOT::Math::Similarity(ca2cu, state.kinematicParametersError().matrix());
0255
0256   //propagation jacobian
0257   AnalyticalCurvilinearJacobian prop(inPar, xPerigee, pPerigee, s);
0258   AlgebraicMatrix66 pr;
0259   pr(5, 5) = 1;
0260   pr.Place_at(prop.jacobian(), 0, 0);
0261
0262   //transportation
0263   AlgebraicSymMatrix66 cov2 = ROOT::Math::Similarity(pr, cov1);
0264
0265   //now geting back to 7-parametrization from curvilinear
0266   cov = ROOT::Math::Similarity(cu2ca, cov2);
0267
0268   // FreeTrajectoryState fts(fPar);
0269
0270   //return parameters as a kiematic state
0271   KinematicParameters resPar(par);
0272   KinematicParametersError resEr(cov);
0273   return KinematicState(resPar, resEr, state.particleCharge(), state.magneticField());
0274
0275   //return  KinematicState(fts,state.mass(), cov(6,6));
0276 }