File indexing completed on 2023-03-17 11:23:24
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
0051 res(1, 1 + i * 7) = 2 * a * sumPy;
0052
0053
0054 res(1, 2 + i * 7) = -2 * a * sumPx;
0055
0056
0057 res(1, 3 + i * 7) = 0.;
0058
0059
0060 res(1, 4 + i * 7) =
0061 2 * sumEnergy / states[i].kinematicParameters().energy() * states[i].kinematicParameters()(3) - 2 * sumPx;
0062
0063
0064 res(1, 5 + i * 7) =
0065 2 * sumEnergy / states[i].kinematicParameters().energy() * states[i].kinematicParameters()(4) - 2 * sumPy;
0066
0067
0068 res(1, 6 + i * 7) =
0069 2 * sumEnergy / states[i].kinematicParameters().energy() * states[i].kinematicParameters()(5) - 2 * sumPz;
0070
0071
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
0096 res(1, 1) = -2 * sumPy * sumA;
0097
0098
0099 res(1, 2) = 2 * sumPx * sumA;
0100
0101
0102 res(1, 3) = 0.;
0103
0104 return res;
0105 }