Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-12 03:09:43

0001 #include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h"
0002 #include "RecoVertex/VertexTools/interface/PerigeeRefittedTrackState.h"
0003 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
0004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0005 #include "MagneticField/Engine/interface/MagneticField.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/Likely.h"
0008 
0009 void PerigeeLinearizedTrackState::computeJacobians() const {
0010   GlobalPoint paramPt(theLinPoint);
0011 
0012   thePredState = builder(theTSOS, paramPt);
0013   if UNLIKELY (!thePredState.isValid())
0014     return;
0015 
0016   double field = theTrack.field()->inInverseGeV(thePredState.theState().position()).z();
0017 
0018   if ((std::abs(theCharge) < 1e-5) || (fabs(field) < 1.e-10)) {
0019     //neutral track
0020     computeNeutralJacobians();
0021   } else {
0022     //charged track
0023     computeChargedJacobians();
0024   }
0025 
0026   jacobiansAvailable = true;
0027 }
0028 
0029 bool PerigeeLinearizedTrackState::operator==(LinearizedTrackState<5>& other) const {
0030   const PerigeeLinearizedTrackState* otherP = dynamic_cast<const PerigeeLinearizedTrackState*>(&other);
0031   if (otherP == nullptr) {
0032     throw VertexException("PerigeeLinearizedTrackState: don't know how to compare myself to non-perigee track state");
0033   }
0034   return (otherP->track() == theTrack);
0035 }
0036 
0037 bool PerigeeLinearizedTrackState::operator==(ReferenceCountingPointer<LinearizedTrackState<5> >& other) const {
0038   const PerigeeLinearizedTrackState* otherP = dynamic_cast<const PerigeeLinearizedTrackState*>(other.get());
0039   if (otherP == nullptr) {
0040     throw VertexException("PerigeeLinearizedTrackState: don't know how to compare myself to non-perigee track state");
0041   }
0042   return (otherP->track() == theTrack);
0043 }
0044 
0045 PerigeeLinearizedTrackState::RefCountedLinearizedTrackState PerigeeLinearizedTrackState::stateWithNewLinearizationPoint(
0046     const GlobalPoint& newLP) const {
0047   return RefCountedLinearizedTrackState(new PerigeeLinearizedTrackState(newLP, track(), theTSOS));
0048 }
0049 
0050 PerigeeLinearizedTrackState::RefCountedRefittedTrackState PerigeeLinearizedTrackState::createRefittedTrackState(
0051     const GlobalPoint& vertexPosition,
0052     const AlgebraicVector3& vectorParameters,
0053     const AlgebraicSymMatrix66& covarianceMatrix) const {
0054   TrajectoryStateClosestToPoint refittedTSCP = PerigeeConversions::trajectoryStateClosestToPoint(
0055       vectorParameters, vertexPosition, charge(), covarianceMatrix, theTrack.field());
0056   return RefCountedRefittedTrackState(new PerigeeRefittedTrackState(refittedTSCP, vectorParameters));
0057 }
0058 
0059 std::vector<PerigeeLinearizedTrackState::RefCountedLinearizedTrackState> PerigeeLinearizedTrackState::components()
0060     const {
0061   std::vector<RefCountedLinearizedTrackState> result;
0062   result.reserve(1);
0063   result.push_back(RefCountedLinearizedTrackState(const_cast<PerigeeLinearizedTrackState*>(this)));
0064   return result;
0065 }
0066 
0067 AlgebraicVector5 PerigeeLinearizedTrackState::refittedParamFromEquation(
0068     const RefCountedRefittedTrackState& theRefittedState) const {
0069   auto p = theRefittedState->position();
0070   AlgebraicVector3 vertexPosition(p.x(), p.y(), p.z());
0071   AlgebraicVector3 momentum = theRefittedState->momentumVector();
0072   if ((momentum(2) * predictedStateMomentumParameters()(2) < 0) && (fabs(momentum(2)) > M_PI / 2)) {
0073     if (predictedStateMomentumParameters()(2) < 0.)
0074       momentum(2) -= 2 * M_PI;
0075     if (predictedStateMomentumParameters()(2) > 0.)
0076       momentum(2) += 2 * M_PI;
0077   }
0078   AlgebraicVectorN param = constantTerm() + positionJacobian() * vertexPosition + momentumJacobian() * momentum;
0079   if (param(2) > M_PI)
0080     param(2) -= 2 * M_PI;
0081   if (param(2) < -M_PI)
0082     param(2) += 2 * M_PI;
0083 
0084   return param;
0085 }
0086 
0087 void PerigeeLinearizedTrackState::checkParameters(AlgebraicVector5& parameters) const {
0088   if (parameters(2) > M_PI)
0089     parameters(2) -= 2 * M_PI;
0090   if (parameters(2) < -M_PI)
0091     parameters(2) += 2 * M_PI;
0092 }
0093 
0094 void PerigeeLinearizedTrackState::computeChargedJacobians() const {
0095   GlobalPoint paramPt(theLinPoint);
0096   //tarjectory parameters
0097   double field = theTrack.field()->inInverseGeV(thePredState.theState().position()).z();
0098   double signTC = -theCharge;
0099 
0100   double thetaAtEP = thePredState.perigeeParameters().theta();
0101   double phiAtEP = thePredState.perigeeParameters().phi();
0102   double ptAtEP = thePredState.pt();
0103   double transverseCurvatureAtEP = field / ptAtEP * signTC;
0104 
0105   double x_v = thePredState.theState().position().x();
0106   double y_v = thePredState.theState().position().y();
0107   double z_v = thePredState.theState().position().z();
0108   double X = x_v - paramPt.x() - sin(phiAtEP) / transverseCurvatureAtEP;
0109   double Y = y_v - paramPt.y() + cos(phiAtEP) / transverseCurvatureAtEP;
0110   double SS = X * X + Y * Y;
0111   double S = sqrt(SS);
0112 
0113   // The track parameters at the expansion point
0114 
0115   theExpandedParams[0] = transverseCurvatureAtEP;
0116   theExpandedParams[1] = thetaAtEP;
0117   theExpandedParams[3] = 1 / transverseCurvatureAtEP - signTC * S;
0118   double phiFEP;
0119   if (std::abs(X) > std::abs(Y)) {
0120     double signX = (X > 0.0 ? +1.0 : -1.0);
0121     phiFEP = -signTC * signX * acos(signTC * Y / S);
0122   } else {
0123     phiFEP = asin(-signTC * X / S);
0124     if ((signTC * Y) < 0.0)
0125       phiFEP = M_PI - phiFEP;
0126   }
0127   if (phiFEP > M_PI)
0128     phiFEP -= 2 * M_PI;
0129   theExpandedParams[2] = phiFEP;
0130   theExpandedParams[4] =
0131       z_v - paramPt.z() - (phiAtEP - theExpandedParams[2]) / tan(thetaAtEP) / transverseCurvatureAtEP;
0132 
0133   // The Jacobian: (all at the expansion point)
0134   // [i,j]
0135   // i = 0: rho , 1: theta, 2: phi_p, 3: epsilon, 4: z_p
0136   // j = 0: x_v, 1: y_v, 2: z_v
0137 
0138   thePositionJacobian(2, 0) = -Y / (SS);
0139   thePositionJacobian(2, 1) = X / (SS);
0140   thePositionJacobian(3, 0) = -signTC * X / S;
0141   thePositionJacobian(3, 1) = -signTC * Y / S;
0142   thePositionJacobian(4, 0) = thePositionJacobian(2, 0) / tan(thetaAtEP) / transverseCurvatureAtEP;
0143   thePositionJacobian(4, 1) = thePositionJacobian(2, 1) / tan(thetaAtEP) / transverseCurvatureAtEP;
0144   thePositionJacobian(4, 2) = 1;
0145 
0146   // [i,j]
0147   // i = 0: rho , 1: theta, 2: phi_p, 3: epsilon, 4: z_p
0148   // j = 0: rho, 1: theta, 2: phi_v
0149   theMomentumJacobian(0, 0) = 1;
0150   theMomentumJacobian(1, 1) = 1;
0151 
0152   theMomentumJacobian(2, 0) =
0153       -(X * cos(phiAtEP) + Y * sin(phiAtEP)) / (SS * transverseCurvatureAtEP * transverseCurvatureAtEP);
0154 
0155   theMomentumJacobian(2, 2) = (Y * cos(phiAtEP) - X * sin(phiAtEP)) / (SS * transverseCurvatureAtEP);
0156 
0157   theMomentumJacobian(3, 0) =
0158       (signTC * (Y * cos(phiAtEP) - X * sin(phiAtEP)) / S - 1) / (transverseCurvatureAtEP * transverseCurvatureAtEP);
0159 
0160   theMomentumJacobian(3, 2) = signTC * (X * cos(phiAtEP) + Y * sin(phiAtEP)) / (S * transverseCurvatureAtEP);
0161 
0162   theMomentumJacobian(4, 0) =
0163       (phiAtEP - theExpandedParams[2]) / tan(thetaAtEP) / (transverseCurvatureAtEP * transverseCurvatureAtEP) +
0164       theMomentumJacobian(2, 0) / tan(thetaAtEP) / transverseCurvatureAtEP;
0165 
0166   theMomentumJacobian(4, 1) =
0167       (phiAtEP - theExpandedParams[2]) * (1 + 1 / (tan(thetaAtEP) * tan(thetaAtEP))) / transverseCurvatureAtEP;
0168 
0169   theMomentumJacobian(4, 2) = (theMomentumJacobian(2, 2) - 1) / tan(thetaAtEP) / transverseCurvatureAtEP;
0170 
0171   // And finally the residuals:
0172 
0173   auto p = thePredState.theState().position();
0174   AlgebraicVector3 expansionPoint(p.x(), p.y(), p.z());
0175   AlgebraicVector3 momentumAtExpansionPoint(transverseCurvatureAtEP, thetaAtEP, phiAtEP);
0176 
0177   theConstantTerm = AlgebraicVector5(theExpandedParams - thePositionJacobian * expansionPoint -
0178                                      theMomentumJacobian * momentumAtExpansionPoint);
0179 }
0180 
0181 void PerigeeLinearizedTrackState::computeNeutralJacobians() const {
0182   GlobalPoint paramPt(theLinPoint);
0183 
0184   //tarjectory parameters
0185   double thetaAtEP = thePredState.theState().momentum().theta();
0186   double phiAtEP = thePredState.theState().momentum().phi();
0187   double ptAtEP = thePredState.theState().momentum().perp();
0188 
0189   double x_v = thePredState.theState().position().x();
0190   double y_v = thePredState.theState().position().y();
0191   double z_v = thePredState.theState().position().z();
0192   double X = x_v - paramPt.x();
0193   double Y = y_v - paramPt.y();
0194 
0195   // The track parameters at the expansion point
0196 
0197   theExpandedParams(0) = 1 / ptAtEP;
0198   theExpandedParams(1) = thetaAtEP;
0199   theExpandedParams(2) = phiAtEP;
0200   theExpandedParams(3) = X * sin(phiAtEP) - Y * cos(phiAtEP);
0201   theExpandedParams(4) = z_v - paramPt.z() - (X * cos(phiAtEP) + Y * sin(phiAtEP)) / tan(thetaAtEP);
0202 
0203   // The Jacobian: (all at the expansion point)
0204   // [i,j]
0205   // i = 0: rho = 1/pt , 1: theta, 2: phi_p, 3: epsilon, 4: z_p
0206   // j = 0: x_v, 1: y_v, 2: z_v
0207 
0208   thePositionJacobian(3, 0) = sin(phiAtEP);
0209   thePositionJacobian(3, 1) = -cos(phiAtEP);
0210   thePositionJacobian(4, 0) = -cos(phiAtEP) / tan(thetaAtEP);
0211   thePositionJacobian(4, 1) = -sin(phiAtEP) / tan(thetaAtEP);
0212   thePositionJacobian(4, 2) = 1;
0213 
0214   // [i,j]
0215   // i = 0: rho = 1/pt , 1: theta, 2: phi_p, 3: epsilon, 4: z_p
0216   // j = 0: rho = 1/pt , 1: theta, 2: phi_v
0217 
0218   theMomentumJacobian(0, 0) = 1;
0219   theMomentumJacobian(1, 1) = 1;
0220   theMomentumJacobian(2, 2) = 1;
0221 
0222   theMomentumJacobian(3, 2) = X * cos(phiAtEP) + Y * sin(phiAtEP);
0223 
0224   theMomentumJacobian(4, 1) = theMomentumJacobian(3, 2) * (1 + 1 / (tan(thetaAtEP) * tan(thetaAtEP)));
0225 
0226   theMomentumJacobian(4, 2) = (X * sin(phiAtEP) - Y * cos(phiAtEP)) / tan(thetaAtEP);
0227 
0228   // And finally the residuals:
0229 
0230   auto p = thePredState.theState().position();
0231   AlgebraicVector3 expansionPoint(p.x(), p.y(), p.z());
0232   AlgebraicVector3 momentumAtExpansionPoint(1. / ptAtEP, thetaAtEP, phiAtEP);
0233 
0234   theConstantTerm = AlgebraicVector5(theExpandedParams - thePositionJacobian * expansionPoint -
0235                                      theMomentumJacobian * momentumAtExpansionPoint);
0236 }