File indexing completed on 2023-10-25 10:03:24
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
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 }
0057
0058 bool TrackKinematicStatePropagator::willPropagateToTheTransversePCA(const KinematicState& state,
0059 const GlobalPoint& point) const {
0060 if (state.particleCharge() == 0.)
0061 return true;
0062
0063
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
0078
0079
0080
0081
0082
0083
0084
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
0120
0121
0122
0123 GlobalTrajectoryParameters fPar(nPosition, nMomentum, state.particleCharge(), state.magneticField());
0124
0125
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
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
0140
0141
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
0155
0156
0157 AlgebraicSymMatrix66 cov1 = ROOT::Math::Similarity(ca2cu, state.kinematicParametersError().matrix());
0158
0159
0160
0161
0162
0163
0164
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
0171 AlgebraicSymMatrix66 cov2 = ROOT::Math::Similarity(pr, cov1);
0172
0173
0174 cov = ROOT::Math::Similarity(cu2ca, cov2);
0175
0176
0177
0178
0179
0180
0181
0182 KinematicParameters resPar(par);
0183 KinematicParametersError resEr(cov);
0184
0185 return KinematicState(resPar, resEr, state.particleCharge(), state.magneticField());
0186
0187
0188
0189
0190
0191
0192 }
0193 }
0194
0195 KinematicState TrackKinematicStatePropagator::propagateToTheTransversePCANeutral(
0196 const KinematicState& state, const GlobalPoint& referencePoint) const {
0197
0198 AlgebraicVector7 par;
0199 AlgebraicSymMatrix77 cov;
0200
0201
0202 GlobalTrajectoryParameters const& inPar = state.trajectoryParameters();
0203
0204
0205
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
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
0233 par(6) = state.mass();
0234
0235
0236
0237
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
0252
0253
0254 AlgebraicSymMatrix66 cov1 = ROOT::Math::Similarity(ca2cu, state.kinematicParametersError().matrix());
0255
0256
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
0263 AlgebraicSymMatrix66 cov2 = ROOT::Math::Similarity(pr, cov1);
0264
0265
0266 cov = ROOT::Math::Similarity(cu2ca, cov2);
0267
0268
0269
0270
0271 KinematicParameters resPar(par);
0272 KinematicParametersError resEr(cov);
0273 return KinematicState(resPar, resEr, state.particleCharge(), state.magneticField());
0274
0275
0276 }