File indexing completed on 2024-04-06 12:29:13
0001 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicPerigeeConversions.h"
0002 #include "TrackingTools/TrajectoryState/interface/PerigeeConversions.h"
0003
0004 ExtendedPerigeeTrajectoryParameters KinematicPerigeeConversions::extendedPerigeeFromKinematicParameters(
0005 const KinematicState& state, const GlobalPoint& point) const {
0006
0007
0008 AlgebraicVector6 res;
0009 res(5) = state.mass();
0010
0011 GlobalVector impactDistance = state.globalPosition() - point;
0012 double theta = state.globalMomentum().theta();
0013 double phi = state.globalMomentum().phi();
0014 double pt = state.globalMomentum().transverse();
0015
0016 double field = state.magneticFieldInInverseGeV().z();
0017
0018
0019
0020
0021
0022 double positiveMomentumPhi = ((phi > 0) ? phi : (2 * M_PI + phi));
0023 double positionPhi = impactDistance.phi();
0024 double positivePositionPhi = ((positionPhi > 0) ? positionPhi : (2 * M_PI + positionPhi));
0025 double phiDiff = positiveMomentumPhi - positivePositionPhi;
0026 if (phiDiff < 0.0)
0027 phiDiff += (2 * M_PI);
0028 double signEpsilon = ((phiDiff > M_PI) ? -1.0 : 1.0);
0029
0030 double epsilon =
0031 signEpsilon * sqrt(impactDistance.x() * impactDistance.x() + impactDistance.y() * impactDistance.y());
0032
0033 double signTC = -state.particleCharge();
0034 bool isCharged = (signTC != 0);
0035
0036 if (isCharged) {
0037 res(0) = field / pt * signTC;
0038 } else {
0039 res(0) = 1 / pt;
0040 }
0041
0042 res(1) = theta;
0043 res(2) = phi;
0044 res(3) = epsilon;
0045 res(4) = impactDistance.z();
0046 return ExtendedPerigeeTrajectoryParameters(res, state.particleCharge());
0047 }
0048
0049 KinematicParameters KinematicPerigeeConversions::kinematicParametersFromExPerigee(
0050 const ExtendedPerigeeTrajectoryParameters& pr, const GlobalPoint& point, const MagneticField* field) const {
0051 AlgebraicVector7 par;
0052 AlgebraicVector6 theVector = pr.vector();
0053 double pt;
0054 if (pr.charge() != 0) {
0055 pt = std::abs(field->inInverseGeV(point).z() / theVector(1));
0056 } else {
0057 pt = 1 / theVector(1);
0058 }
0059 par(6) = theVector[5];
0060 par(0) = theVector[3] * sin(theVector[2]) + point.x();
0061 par(1) = -theVector[3] * cos(theVector[2]) + point.y();
0062 par(2) = theVector[4] + point.z();
0063 par(3) = cos(theVector[2]) * pt;
0064 par(4) = sin(theVector[2]) * pt;
0065 par(5) = pt / tan(theVector[1]);
0066
0067 return KinematicParameters(par);
0068 }
0069
0070 KinematicState KinematicPerigeeConversions::kinematicState(const AlgebraicVector4& momentum,
0071 const GlobalPoint& referencePoint,
0072 const TrackCharge& charge,
0073 const AlgebraicSymMatrix77& theCovarianceMatrix,
0074 const MagneticField* field) const {
0075 AlgebraicMatrix77 param2cart = jacobianParameters2Kinematic(momentum, referencePoint, charge, field);
0076 AlgebraicSymMatrix77 kinematicErrorMatrix = ROOT::Math::Similarity(param2cart, theCovarianceMatrix);
0077
0078
0079 KinematicParametersError kinematicParamError(kinematicErrorMatrix);
0080 AlgebraicVector7 par;
0081 AlgebraicVector4 mm = momentumFromPerigee(momentum, referencePoint, charge, field);
0082 par(0) = referencePoint.x();
0083 par(1) = referencePoint.y();
0084 par(2) = referencePoint.z();
0085 par(3) = mm(0);
0086 par(4) = mm(1);
0087 par(5) = mm(2);
0088 par(6) = mm(3);
0089 KinematicParameters kPar(par);
0090 return KinematicState(kPar, kinematicParamError, charge, field);
0091 }
0092
0093 AlgebraicMatrix77 KinematicPerigeeConversions::jacobianParameters2Kinematic(const AlgebraicVector4& momentum,
0094 const GlobalPoint& referencePoint,
0095 const TrackCharge& charge,
0096 const MagneticField* field) const {
0097 AlgebraicMatrix66 param2cart = PerigeeConversions::jacobianParameters2Cartesian(
0098 AlgebraicVector3(momentum[0], momentum[1], momentum[2]), referencePoint, charge, field);
0099 AlgebraicMatrix77 frameTransJ;
0100 for (int i = 0; i < 6; ++i)
0101 for (int j = 0; j < 6; ++j)
0102 frameTransJ(i, j) = param2cart(i, j);
0103 frameTransJ(6, 6) = 1;
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123 return frameTransJ;
0124 }
0125
0126
0127 AlgebraicVector4 KinematicPerigeeConversions::momentumFromPerigee(const AlgebraicVector4& momentum,
0128 const GlobalPoint& referencePoint,
0129 const TrackCharge& ch,
0130 const MagneticField* field) const {
0131 AlgebraicVector4 mm;
0132 double pt;
0133 if (ch != 0) {
0134
0135 pt = std::abs(field->inInverseGeV(referencePoint).z() / momentum[0]);
0136 } else {
0137 pt = 1 / momentum[0];
0138 }
0139 mm(0) = cos(momentum[2]) * pt;
0140 mm(1) = sin(momentum[2]) * pt;
0141 mm(2) = pt / tan(momentum[1]);
0142 mm(3) = momentum[3];
0143 return mm;
0144 }