KinematicFit/src/KinematicConstrainedVertexFitter.cc

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   //sorting out the input particles
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   // linearization point
0061   // (only compute it using the linearization point finder if no point was passed to the fit function):
0062   GlobalPoint linPoint;
0063   if (pt != nullptr) {
0064     linPoint = *pt;
0065   } else {
0066     linPoint = finder->getLinearizationPoint(fStates);
0067   }
0068
0069   //initial parameters:
0070   int vSize = particles.size();
0071   AlgebraicVector inPar(3 + 7 * vSize, 0);
0072
0073   //final parameters
0074   AlgebraicVector finPar(3 + 7 * vSize, 0);
0075
0076   //initial covariance
0077   AlgebraicMatrix inCov(3 + 7 * vSize, 3 + 7 * vSize, 0);
0078
0079   //making initial vector of parameters and initial particle-related covariance
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   //initial vertex error matrix components (huge error method)
0099   //and vertex related initial vector components
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   //constraint equations value and number of iterations
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   //iterarions over the updator: each time updated parameters
0123   //are taken as new linearization point
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   // std::cout << "old full cov matrix" << std::endl;
0185   // std::cout << refCCov << std::endl;
0186
0187   // cout<<"number of relinearizations "<<nit<<endl;
0188   // cout<<"value obtained: "<<eq<<endl;
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; }