KinematicFitPrimitives/src/KinematicPerigeeConversions.cc

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   //making an extended perigee parametrization
0007   //out of kinematic state and point defined
0008   AlgebraicVector6 res;
0009   res(5) = state.mass();
0010   //  AlgebraicVector7 par = state.kinematicParameters().vector();
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   //  double field = MagneticField::inGeVPerCentimeter(state.globalPosition()).z();
0016   double field = state.magneticFieldInInverseGeV().z();
0017   // double signTC = -state.particleCharge();
0018   // double transverseCurvature = field/pt*signTC;
0019
0020   //making a proper sign for epsilon
0021   // FIXME use scalar product...
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   //  kinematicErrorMatrix.assign(param2cart*theCovarianceMatrix*param2cart.T());
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   //   double factor = 1;
0106   //   if (charge != 0){
0107   //    double field = TrackingTools::FakeField::Field::inGeVPerCentimeter(referencePoint).z();
0108   //    factor =  -field*charge;
0109   //    }
0110   //   AlgebraicMatrix frameTransJ(7, 7, 0);
0111   //   frameTransJ[0][0] = 1;
0112   //   frameTransJ[1][1] = 1;
0113   //   frameTransJ[2][2] = 1;
0114   //   frameTransJ[6][6] = 1;
0115   //   frameTransJ[3][3] = - factor * cos(momentum[2])  / (momentum[0]*momentum[0]);
0116   //   frameTransJ[4][3] = - factor * sin(momentum[2])  / (momentum[0]*momentum[0]);
0117   //   frameTransJ[5][3] = - factor / tan(momentum[1])  / (momentum[0]*momentum[0]);
0118   //
0119   //   frameTransJ[3][5] = - factor * sin(momentum[1]) / (momentum[0]);
0120   //   frameTransJ[4][5] =   factor * cos(momentum[1])  / (momentum[0]);
0121   //   frameTransJ[5][4] = -factor/ (momentum[0]*sin(momentum[1])*sin(momentum[1]));
0122
0123   return frameTransJ;
0124 }
0125
0126 // Cartesian (px,py,px,m) from extended perigee
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     //   pt = abs(MagneticField::inGeVPerCentimeter(referencePoint).z() / momentum[0]);
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 }