KinematicFit/src/MultiTrackPointingKinematicConstraint.cc

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