File indexing completed on 2024-04-06 12:29:11
0001 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
0002 #include "RecoVertex/KinematicFit/interface/InputSort.h"
0003 #include "RecoVertex/LinearizationPointFinders/interface/DefaultLinearizationPointFinder.h"
0004 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicVertexFactory.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter() {
0008 finder = new DefaultLinearizationPointFinder();
0009 vCons = new VertexKinematicConstraint();
0010 updator = new KinematicConstrainedVertexUpdator();
0011 tBuilder = new ConstrainedTreeBuilder;
0012 defaultParameters();
0013 iterations = -1;
0014 csum = -1000.0;
0015 }
0016
0017 KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter(const LinearizationPointFinder &fnd) {
0018 finder = fnd.clone();
0019 vCons = new VertexKinematicConstraint();
0020 updator = new KinematicConstrainedVertexUpdator();
0021 tBuilder = new ConstrainedTreeBuilder;
0022 defaultParameters();
0023 iterations = -1;
0024 csum = -1000.0;
0025 }
0026
0027 KinematicConstrainedVertexFitter::~KinematicConstrainedVertexFitter() {
0028 delete finder;
0029 delete vCons;
0030 delete updator;
0031 delete tBuilder;
0032 }
0033
0034 void KinematicConstrainedVertexFitter::setParameters(const edm::ParameterSet &pSet) {
0035 theMaxDelta = pSet.getParameter<double>("maxDelta");
0036 theMaxStep = pSet.getParameter<int>("maxNbrOfIterations");
0037 theMaxReducedChiSq = pSet.getParameter<double>("maxReducedChiSq");
0038 theMinChiSqImprovement = pSet.getParameter<double>("minChiSqImprovement");
0039 }
0040
0041 void KinematicConstrainedVertexFitter::defaultParameters() {
0042 theMaxDelta = 0.01;
0043 theMaxStep = 1000;
0044 theMaxReducedChiSq = 225.;
0045 theMinChiSqImprovement = 50.;
0046 }
0047
0048 RefCountedKinematicTree KinematicConstrainedVertexFitter::fit(const std::vector<RefCountedKinematicParticle> &part,
0049 MultiTrackKinematicConstraint *cs,
0050 GlobalPoint *pt) {
0051 if (part.size() < 2)
0052 throw VertexException("KinematicConstrainedVertexFitter::input states are less than 2");
0053
0054
0055 InputSort iSort;
0056 std::pair<std::vector<RefCountedKinematicParticle>, std::vector<FreeTrajectoryState> > input = iSort.sort(part);
0057 const std::vector<RefCountedKinematicParticle> &particles = input.first;
0058 const std::vector<FreeTrajectoryState> &fStates = input.second;
0059
0060
0061
0062 GlobalPoint linPoint;
0063 if (pt != nullptr) {
0064 linPoint = *pt;
0065 } else {
0066 linPoint = finder->getLinearizationPoint(fStates);
0067 }
0068
0069
0070 int vSize = particles.size();
0071 AlgebraicVector inPar(3 + 7 * vSize, 0);
0072
0073
0074 AlgebraicVector finPar(3 + 7 * vSize, 0);
0075
0076
0077 AlgebraicMatrix inCov(3 + 7 * vSize, 3 + 7 * vSize, 0);
0078
0079
0080 int nSt = 0;
0081 std::vector<KinematicState> inStates;
0082 for (std::vector<RefCountedKinematicParticle>::const_iterator i = particles.begin(); i != particles.end(); i++) {
0083 KinematicState state = (*i)->stateAtPoint(linPoint);
0084 if (!state.isValid()) {
0085 LogDebug("KinematicConstrainedVertexFitter") << "State is invalid at point: " << linPoint << std::endl;
0086 return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0087 }
0088 AlgebraicVector prPar = asHepVector<7>(state.kinematicParameters().vector());
0089 for (int j = 1; j < 8; j++) {
0090 inPar(3 + 7 * nSt + j) = prPar(j);
0091 }
0092 AlgebraicSymMatrix l_cov = asHepMatrix<7>(state.kinematicParametersError().matrix());
0093 inCov.sub(4 + 7 * nSt, 4 + 7 * nSt, l_cov);
0094 inStates.push_back(state);
0095 ++nSt;
0096 }
0097
0098
0099
0100 double in_er = 100.;
0101 inCov(1, 1) = in_er;
0102 inCov(2, 2) = in_er;
0103 inCov(3, 3) = in_er;
0104
0105 inPar(1) = linPoint.x();
0106 inPar(2) = linPoint.y();
0107 inPar(3) = linPoint.z();
0108
0109
0110 double eq;
0111 int nit = 0;
0112 iterations = 0;
0113 csum = 0.0;
0114
0115 std::vector<KinematicState> lStates = inStates;
0116 GlobalPoint lPoint = linPoint;
0117 RefCountedKinematicVertex rVtx;
0118 AlgebraicMatrix refCCov;
0119
0120 double chisq = 1e6;
0121 bool convergence = false;
0122
0123
0124 do {
0125 eq = 0.;
0126 std::pair<std::pair<std::vector<KinematicState>, AlgebraicMatrix>, RefCountedKinematicVertex> lRes =
0127 updator->update(inPar, inCov, lStates, lPoint, cs);
0128
0129 const std::vector<KinematicState> &newStates = lRes.first.first;
0130
0131 if (particles.size() != newStates.size()) {
0132 LogDebug("KinematicConstrainedVertexFitter") << "updator failure\n";
0133 return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0134 }
0135
0136 rVtx = lRes.second;
0137
0138 double newchisq = rVtx->chiSquared();
0139 if (nit > 2 && newchisq > theMaxReducedChiSq * rVtx->degreesOfFreedom() &&
0140 (newchisq - chisq) > (-theMinChiSqImprovement)) {
0141 LogDebug("KinematicConstrainedVertexFitter") << "bad chisq and insufficient improvement, bailing\n";
0142 return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0143 }
0144 chisq = newchisq;
0145
0146 const GlobalPoint &newPoint = rVtx->position();
0147
0148 double maxDelta = 0.0;
0149
0150 double deltapos[3];
0151 deltapos[0] = newPoint.x() - lPoint.x();
0152 deltapos[1] = newPoint.y() - lPoint.y();
0153 deltapos[2] = newPoint.z() - lPoint.z();
0154 for (int i = 0; i < 3; ++i) {
0155 double delta = deltapos[i] * deltapos[i] / rVtx->error().matrix()(i, i);
0156 if (delta > maxDelta)
0157 maxDelta = delta;
0158 }
0159
0160 for (std::vector<KinematicState>::const_iterator itold = lStates.begin(), itnew = newStates.begin();
0161 itnew != newStates.end();
0162 ++itold, ++itnew) {
0163 for (int i = 0; i < 7; ++i) {
0164 double deltapar = itnew->kinematicParameters()(i) - itold->kinematicParameters()(i);
0165 double delta = deltapar * deltapar / itnew->kinematicParametersError().matrix()(i, i);
0166 if (delta > maxDelta)
0167 maxDelta = delta;
0168 }
0169 }
0170
0171 lStates = newStates;
0172 lPoint = newPoint;
0173
0174 refCCov = lRes.first.second;
0175 nit++;
0176 convergence = maxDelta < theMaxDelta || (nit == theMaxStep && maxDelta < 4.0 * theMaxDelta);
0177
0178 } while (nit < theMaxStep && !convergence);
0179
0180 if (!convergence) {
0181 return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0182 }
0183
0184
0185
0186
0187
0188
0189 iterations = nit;
0190 csum = eq;
0191
0192 return tBuilder->buildTree(particles, lStates, rVtx, refCCov);
0193 }
0194
0195 int KinematicConstrainedVertexFitter::getNit() const { return iterations; }
0196
0197 float KinematicConstrainedVertexFitter::getCSum() const { return csum; }