KinematicFit/src/MultiTrackMassKinematicConstraint.cc

0001 #include "RecoVertex/KinematicFit/interface/MultiTrackMassKinematicConstraint.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003
0004 AlgebraicVector MultiTrackMassKinematicConstraint::value(const std::vector<KinematicState>& states,
0005                                                          const GlobalPoint& point) const {
0006   if (states.size() < nPart)
0007     throw VertexException("MultiTrackMassKinematicConstraint::not enough states given");
0008
0009   double sumEnergy = 0, sumPx = 0, sumPy = 0., sumPz = 0.;
0010
0011   double a;
0012   for (unsigned int i = 0; i < nPart; ++i) {
0013     a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();
0014
0015     sumEnergy += states[i].kinematicParameters().energy();
0016     sumPx += states[i].kinematicParameters()(3) - a * (point.y() - states[i].kinematicParameters()(1));
0017     sumPy += states[i].kinematicParameters()(4) + a * (point.x() - states[i].kinematicParameters()(0));
0018     sumPz += states[i].kinematicParameters()(5);
0019   }
0020
0021   double j_m = sumPx * sumPx + sumPy * sumPy + sumPz * sumPz;
0022
0023   AlgebraicVector res(1, 0);
0024   res(1) = sumEnergy * sumEnergy - j_m - mass * mass;
0025   return res;
0026 }
0027
0028 AlgebraicMatrix MultiTrackMassKinematicConstraint::parametersDerivative(const std::vector<KinematicState>& states,
0029                                                                         const GlobalPoint& point) const {
0030   if (states.size() < nPart)
0031     throw VertexException("MultiTrackMassKinematicConstraint::not enough states given");
0032
0033   AlgebraicMatrix res(1, states.size() * 7, 0);
0034
0035   double sumEnergy = 0, sumPx = 0, sumPy = 0., sumPz = 0.;
0036
0037   double a;
0038   for (unsigned int i = 0; i < nPart; ++i) {
0039     a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();
0040
0041     sumEnergy += states[i].kinematicParameters().energy();
0042     sumPx += states[i].kinematicParameters()(3) - a * (point.y() - states[i].kinematicParameters()(1));
0043     sumPy += states[i].kinematicParameters()(4) + a * (point.x() - states[i].kinematicParameters()(0));
0044     sumPz += states[i].kinematicParameters()(5);
0045   }
0046
0047   for (unsigned int i = 0; i < nPart; ++i) {
0048     a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();
0049
0050     //x derivatives:
0051     res(1, 1 + i * 7) = 2 * a * sumPy;
0052
0053     //y derivatives:
0054     res(1, 2 + i * 7) = -2 * a * sumPx;
0055
0056     //z components:
0057     res(1, 3 + i * 7) = 0.;
0058
0059     //px components:
0060     res(1, 4 + i * 7) =
0061         2 * sumEnergy / states[i].kinematicParameters().energy() * states[i].kinematicParameters()(3) - 2 * sumPx;
0062
0063     //py components:
0064     res(1, 5 + i * 7) =
0065         2 * sumEnergy / states[i].kinematicParameters().energy() * states[i].kinematicParameters()(4) - 2 * sumPy;
0066
0067     //pz1 components:
0068     res(1, 6 + i * 7) =
0069         2 * sumEnergy / states[i].kinematicParameters().energy() * states[i].kinematicParameters()(5) - 2 * sumPz;
0070
0071     //mass components:
0072     res(1, 7 + i * 7) =
0073         2 * states[i].kinematicParameters().mass() * sumEnergy / states[i].kinematicParameters().energy();
0074   }
0075   return res;
0076 }
0077
0078 AlgebraicMatrix MultiTrackMassKinematicConstraint::positionDerivative(const std::vector<KinematicState>& states,
0079                                                                       const GlobalPoint& point) const {
0080   AlgebraicMatrix res(1, 3, 0);
0081   if (states.size() < nPart)
0082     throw VertexException("MultiTrackMassKinematicConstraint::not enough states given");
0083
0084   double sumA = 0, sumPx = 0, sumPy = 0.;
0085
0086   double a;
0087   for (unsigned int i = 0; i < nPart; ++i) {
0088     a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();
0089     sumA += a;
0090
0091     sumPx += states[i].kinematicParameters()(3) - a * (point.y() - states[i].kinematicParameters()(1));
0092     sumPy += states[i].kinematicParameters()(4) + a * (point.x() - states[i].kinematicParameters()(0));
0093   }
0094
0095   //xv component
0096   res(1, 1) = -2 * sumPy * sumA;
0097
0098   //yv component
0099   res(1, 2) = 2 * sumPx * sumA;
0100
0101   //zv component
0102   res(1, 3) = 0.;
0103
0104   return res;
0105 }