Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:17

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