KinematicFit/src/MultiTrackVertexLinkKinematicConstraint.cc

0001 #include "RecoVertex/KinematicFit/interface/MultiTrackVertexLinkKinematicConstraint.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003
0004 AlgebraicVector MultiTrackVertexLinkKinematicConstraint::value(const std::vector<KinematicState>& states,
0005                                                                const GlobalPoint& point) const {
0006   int num = states.size();
0007   if (num < 2)
0008     throw VertexException("MultiTrackVertexLinkKinematicConstraint::value <2 states passed");
0009
0010   //2 equations (for all tracks)
0011   AlgebraicVector vl(2, 0);
0012   double dx = point.x() - refPoint.x();
0013   double dy = point.y() - refPoint.y();
0014   double dz = point.z() - refPoint.z();
0015   double dT = sqrt(pow(dx, 2) + pow(dy, 2));
0016   double ds = sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2));
0017
0018   double pxSum = 0, pySum = 0, pzSum = 0;
0019   double aSum = 0;
0020   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0021     double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
0022     aSum += a;
0023
0024     pxSum += i->kinematicParameters()(3) - a * (point.y() - i->kinematicParameters()(1));
0025     pySum += i->kinematicParameters()(4) + a * (point.x() - i->kinematicParameters()(0));
0026     pzSum += i->kinematicParameters()(5);
0027   }
0028
0029   double pT = sqrt(pow(pxSum, 2) + pow(pySum, 2));
0030   double pSum = sqrt(pow(pxSum, 2) + pow(pySum, 2) + pow(pzSum, 2));
0031
0032   vl(1) = (dT - dx) / dy + (-2 * pT + sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) / (aSum * dT) +
0033           (-pT + pxSum) / pySum;
0034   vl(2) = (ds - dT) / dz + (pT - pSum) / pzSum;
0035
0036   return vl;
0037 }
0038
0039 AlgebraicMatrix MultiTrackVertexLinkKinematicConstraint::parametersDerivative(const std::vector<KinematicState>& states,
0040                                                                               const GlobalPoint& point) const {
0041   int num = states.size();
0042   if (num < 2)
0043     throw VertexException("MultiTrackVertexLinkKinematicConstraint::parametersDerivative <2 states passed");
0044
0045   //2 equations (for all tracks)
0046   AlgebraicMatrix matrix(2, num * 7, 0);  //AlgebraicMatrix starts from 1
0047   double dx = point.x() - refPoint.x();
0048   double dy = point.y() - refPoint.y();
0049   double dT = sqrt(pow(dx, 2) + pow(dy, 2));
0050
0051   double pxSum = 0, pySum = 0, pzSum = 0;
0052   double aSum = 0;
0053   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0054     double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
0055     aSum += a;
0056
0057     pxSum += i->kinematicParameters()(3) - a * (point.y() - i->kinematicParameters()(1));
0058     pySum += i->kinematicParameters()(4) + a * (point.x() - i->kinematicParameters()(0));
0059     pzSum += i->kinematicParameters()(5);
0060   }
0061
0062   double pT = sqrt(pow(pxSum, 2) + pow(pySum, 2));
0063   double pSum = sqrt(pow(pxSum, 2) + pow(pySum, 2) + pow(pzSum, 2));
0064
0065   int col = 0;
0066   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0067     double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
0068
0069     matrix(1, 1 + col * 7) = a * (-(pT / pow(pySum, 2)) + pxSum / pow(pySum, 2) -
0070                                   (4 * pySum) / (aSum * dT * sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) +
0071                                   (1 + (2 * pySum) / (aSum * dT)) / pT);  //dH/dx
0072     matrix(1, 2 + col * 7) =
0073         (a * (aSum * dT * (pT - pxSum) +
0074               2 * (-1 + (2 * pT) / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) * pxSum * pySum)) /
0075         (aSum * dT * pT * pySum);  //dH/dy
0076     //dH/dz=0
0077     matrix(1, 4 + col * 7) =
0078         (aSum * dT * (pT - pxSum) +
0079          2 * (-1 + (2 * pT) / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) * pxSum * pySum) /
0080         (aSum * dT * pT * pySum);  //dH/dpx
0081     matrix(1, 5 + col * 7) = pT / pow(pySum, 2) - pxSum / pow(pySum, 2) +
0082                              (4 * pySum) / (aSum * dT * sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) +
0083                              (-1 - (2 * pySum) / (aSum * dT)) / pT;  //dH/dpy
0084     //dH/dpz=0
0085     //dH/dm=0
0086     matrix(2, 1 + col * 7) = (a * (-pSum + pT) * pySum) / (pSum * pT * pzSum);  //dH/dx
0087     matrix(2, 2 + col * 7) = (a * (pSum - pT) * pxSum) / (pSum * pT * pzSum);   //dH/dy
0088     //dH/dz
0089     matrix(2, 4 + col * 7) = ((-(1 / pSum) + 1 / pT) * pxSum) / pzSum;   //dH/dpx
0090     matrix(2, 5 + col * 7) = ((-(1 / pSum) + 1 / pT) * pySum) / pzSum;   //dH/dpy
0091     matrix(2, 6 + col * 7) = -(1 / pSum) + (pSum - pT) / pow(pzSum, 2);  //dH/dpz
0092     //dH/dm=0
0093
0094     col++;
0095   }
0096
0097   return matrix;
0098 }
0099
0100 AlgebraicMatrix MultiTrackVertexLinkKinematicConstraint::positionDerivative(const std::vector<KinematicState>& states,
0101                                                                             const GlobalPoint& point) const {
0102   int num = states.size();
0103   if (num < 2)
0104     throw VertexException("MultiTrackVertexLinkKinematicConstraint::positionDerivative <2 states passed");
0105
0106   //2 equations (for all tracks)
0107   AlgebraicMatrix matrix(2, 3, 0);
0108   double dx = point.x() - refPoint.x();
0109   double dy = point.y() - refPoint.y();
0110   double dz = point.z() - refPoint.z();
0111   double dT = sqrt(pow(dx, 2) + pow(dy, 2));
0112   double ds = sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2));
0113
0114   double pxSum = 0, pySum = 0, pzSum = 0, aSum = 0;
0115   for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
0116     double a = -i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
0117     aSum += a;
0118
0119     pxSum += i->kinematicParameters()(3) - a * (point.y() - i->kinematicParameters()(1));
0120     pySum += i->kinematicParameters()(4) + a * (point.x() - i->kinematicParameters()(0));
0121     pzSum += i->kinematicParameters()(5);
0122   }
0123   double pT = sqrt(pow(pxSum, 2) + pow(pySum, 2));
0124   double pSum = sqrt(pow(pxSum, 2) + pow(pySum, 2) + pow(pzSum, 2));
0125
0126   matrix(1, 1) =
0127       (-1 + dx / dT) / dy +
0128       (2 * dx * pT * (1 - (2 * pT) / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2)))) / (aSum * pow(dT, 3)) +
0129       aSum * (-(1 / pT) + pT / pow(pySum, 2) - pxSum / pow(pySum, 2)) +
0130       (2 * (-(1 / pT) + 2 / sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) * pySum) / dT;  //dH/dxv
0131   matrix(1, 2) = 1 / dT + (-dT + dx) / pow(dy, 2) -
0132                  (dy * (-2 * pT + sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2)))) / (aSum * pow(dT, 3)) -
0133                  ((-2 + sqrt(4 - (pow(aSum, 2) * pow(dT, 2)) / pow(pT, 2))) * pxSum) / (dT * pT) -
0134                  (aSum * (dy * pow(pT, 2) + aSum * pow(dT, 2) * pxSum)) /
0135                      (dT * pow(pT, 2) * sqrt(-(pow(aSum, 2) * pow(dT, 2)) + 4 * pow(pT, 2))) +
0136                  (aSum * (-pT + pxSum)) / (pT * pySum);  //dH/dyv
0137   //dH/dzv=0
0138   matrix(2, 1) = ((1 / ds - 1 / dT) * dx) / dz + (aSum * (pSum - pT) * pySum) / (pSum * pT * pzSum);  //dH/dxv
0139   matrix(2, 2) = ((1 / ds - 1 / dT) * dy) / dz - (aSum * (pSum - pT) * pxSum) / (pSum * pT * pzSum);  //dH/dyv
0140   matrix(2, 3) = 1 / ds + (-ds + dT) / pow(dz, 2);                                                    //dH/dzv
0141
0142   return matrix;
0143 }
0144
0145 int MultiTrackVertexLinkKinematicConstraint::numberOfEquations() const { return 2; }