KinematicFit/interface/KinematicConstrainedVertexFitterT.h

0001 #ifndef KinematicConstrainedVertexFitterT_H
0002 #define KinematicConstrainedVertexFitterT_H
0003
0004 #include "RecoVertex/KinematicFitPrimitives/interface/RefCountedKinematicTree.h"
0005 #include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraintT.h"
0006 #include "RecoVertex/VertexTools/interface/LinearizationPointFinder.h"
0007 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexUpdatorT.h"
0008 #include "RecoVertex/KinematicFit/interface/VertexKinematicConstraintT.h"
0009
0010 #include "RecoVertex/KinematicFit/interface/ConstrainedTreeBuilderT.h"
0011
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014
0015 /**
0016  * Class fitting the veretx out of set of tracks via usual LMS
0017  * with Lagrange multipliers.
0018  * Additional constraints can be applyed to the tracks during the vertex fit
0019  * (solves non-factorizabele cases). Since the vertex constraint is included by default, do not add a separate
0020  * VertexKinematicConstraint!
0021  * Example: Vertex fit with collinear tracks..
0022  */
0023
0024 template <int nTrk, int nConstraint>
0025 class KinematicConstrainedVertexFitterT {
0026 public:
0027   /**
0028    * Default constructor using LMSLinearizationPointFinder
0029    */
0030   explicit KinematicConstrainedVertexFitterT(const MagneticField* ifield);
0031
0032   /**
0033    * Constructor with user-provided LinearizationPointFinder
0034    */
0035   KinematicConstrainedVertexFitterT(const MagneticField* ifield, const LinearizationPointFinder& fnd);
0036
0037   ~KinematicConstrainedVertexFitterT();
0038
0039   /**
0040    * Configuration through PSet: number of iterations(maxDistance) and
0041    * stopping condition (maxNbrOfIterations)
0042    */
0043
0044   void setParameters(const edm::ParameterSet& pSet);
0045
0046   /**
0047    * Without additional constraint, this will perform a simple
0048    * vertex fit using LMS with Lagrange multipliers method (by definition valid only if nConstraint=0)
0049    */
0050   RefCountedKinematicTree fit(const std::vector<RefCountedKinematicParticle>& part) { return fit(part, 0, 0); }
0051
0052   /**
0053    * LMS with Lagrange multipliers fit of vertex constraint and user-specified constraint.
0054    */
0055   RefCountedKinematicTree fit(const std::vector<RefCountedKinematicParticle>& part,
0056                               MultiTrackKinematicConstraintT<nTrk, nConstraint>* cs) {
0057     return fit(part, cs, nullptr);
0058   };
0059
0060   /**
0061    * LMS with Lagrange multipliers fit of vertex constraint, user-specified constraint and user-specified starting point.
0062    */
0063   RefCountedKinematicTree fit(const std::vector<RefCountedKinematicParticle>& part,
0064                               MultiTrackKinematicConstraintT<nTrk, nConstraint>* cs,
0065                               GlobalPoint* pt);
0066
0067   //return the number of iterations
0068   int getNit() const;
0069   //return the value of the constraint equation
0070   float getCSum() const;
0071
0072 private:
0073   void defaultParameters();
0074
0075   const MagneticField* field;
0076   LinearizationPointFinder* finder;
0077   KinematicConstrainedVertexUpdatorT<nTrk, nConstraint>* updator;
0078   VertexKinematicConstraintT* vCons;
0079   ConstrainedTreeBuilderT* tBuilder;
0080
0081   float theMaxDelta;  //maximum (delta parameter)^2/(sigma parameter)^2 per iteration for convergence
0082   int theMaxStep;
0083   float theMaxReducedChiSq;  //max of initial (after 2 iterations) chisq/dof value
0084   float theMinChiSqImprovement;  //minimum required improvement in chisq to avoid fit termination for cases exceeding theMaxReducedChiSq
0085
0086   int iterations;
0087   float csum;
0088 };
0089
0090 #include "RecoVertex/KinematicFit/interface/InputSort.h"
0091 #include "RecoVertex/LinearizationPointFinders/interface/DefaultLinearizationPointFinder.h"
0092 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicVertexFactory.h"
0093 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0094
0095 template <int nTrk, int nConstraint>
0096 KinematicConstrainedVertexFitterT<nTrk, nConstraint>::KinematicConstrainedVertexFitterT(const MagneticField* ifield)
0097     : field(ifield) {
0098   finder = new DefaultLinearizationPointFinder();
0099   vCons = new VertexKinematicConstraintT();
0100   updator = new KinematicConstrainedVertexUpdatorT<nTrk, nConstraint>();
0101   tBuilder = new ConstrainedTreeBuilderT;
0102   defaultParameters();
0103   iterations = -1;
0104   csum = -1000.0;
0105 }
0106
0107 template <int nTrk, int nConstraint>
0108 KinematicConstrainedVertexFitterT<nTrk, nConstraint>::KinematicConstrainedVertexFitterT(
0109     const MagneticField* ifield, const LinearizationPointFinder& fnd)
0110     : field(ifield) {
0111   finder = fnd.clone();
0112   vCons = new VertexKinematicConstraintT();
0113   updator = new KinematicConstrainedVertexUpdatorT<nTrk, nConstraint>();
0114   tBuilder = new ConstrainedTreeBuilderT;
0115   defaultParameters();
0116   iterations = -1;
0117   csum = -1000.0;
0118 }
0119
0120 template <int nTrk, int nConstraint>
0121 KinematicConstrainedVertexFitterT<nTrk, nConstraint>::~KinematicConstrainedVertexFitterT() {
0122   delete finder;
0123   delete vCons;
0124   delete updator;
0125   delete tBuilder;
0126 }
0127
0128 template <int nTrk, int nConstraint>
0129 void KinematicConstrainedVertexFitterT<nTrk, nConstraint>::setParameters(const edm::ParameterSet& pSet) {
0130   theMaxDelta = pSet.getParameter<double>("maxDelta");
0131   theMaxStep = pSet.getParameter<int>("maxNbrOfIterations");
0132   theMaxReducedChiSq = pSet.getParameter<double>("maxReducedChiSq");
0133   theMinChiSqImprovement = pSet.getParameter<double>("minChiSqImprovement");
0134 }
0135
0136 template <int nTrk, int nConstraint>
0137 void KinematicConstrainedVertexFitterT<nTrk, nConstraint>::defaultParameters() {
0138   theMaxDelta = 0.01;
0139   theMaxStep = 1000;
0140   theMaxReducedChiSq = 225.;
0141   theMinChiSqImprovement = 50.;
0142 }
0143
0144 template <int nTrk, int nConstraint>
0145 RefCountedKinematicTree KinematicConstrainedVertexFitterT<nTrk, nConstraint>::fit(
0146     const std::vector<RefCountedKinematicParticle>& part,
0147     MultiTrackKinematicConstraintT<nTrk, nConstraint>* cs,
0148     GlobalPoint* pt) {
0149   assert(nConstraint == 0 || cs != nullptr);
0150   if (part.size() != nTrk)
0151     throw VertexException("KinematicConstrainedVertexFitterT::input states are not nTrk");
0152
0153   //sorting out the input particles
0154   InputSort iSort;
0155   std::pair<std::vector<RefCountedKinematicParticle>, std::vector<FreeTrajectoryState> > input = iSort.sort(part);
0156   const std::vector<RefCountedKinematicParticle>& particles = input.first;
0157   const std::vector<FreeTrajectoryState>& fStates = input.second;
0158
0159   // linearization point:
0160   GlobalPoint linPoint = (pt != nullptr) ? *pt : finder->getLinearizationPoint(fStates);
0161
0162   //initial parameters:
0163   ROOT::Math::SVector<double, 3 + 7 * nTrk> inPar;   //3+ 7*ntracks
0164   ROOT::Math::SVector<double, 3 + 7 * nTrk> finPar;  //3+ 7*ntracks
0165
0166   ROOT::Math::SMatrix<double, 3 + 7 * nTrk, 3 + 7 * nTrk, ROOT::Math::MatRepSym<double, 3 + 7 * nTrk> > inCov;
0167
0168   //making initial vector of parameters and initial particle-related covariance
0169   int nSt = 0;
0170   std::vector<KinematicState> lStates(nTrk);
0171   for (std::vector<RefCountedKinematicParticle>::const_iterator i = particles.begin(); i != particles.end(); i++) {
0172     lStates[nSt] = (*i)->stateAtPoint(linPoint);
0173     KinematicState const& state = lStates[nSt];
0174     if (!state.isValid()) {
0175       LogDebug("KinematicConstrainedVertexFitter") << "State is invalid at point: " << linPoint << std::endl;
0176       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0177     }
0178     inPar.Place_at(state.kinematicParameters().vector(), 3 + 7 * nSt);
0179     inCov.Place_at(state.kinematicParametersError().matrix(), 3 + 7 * nSt, 3 + 7 * nSt);
0180     ++nSt;
0181   }
0182
0183   //initial vertex error matrix components (huge error method)
0184   //and vertex related initial vector components
0185   double in_er = 100.;
0186   inCov(0, 0) = in_er;
0187   inCov(1, 1) = in_er;
0188   inCov(2, 2) = in_er;
0189
0190   inPar(0) = linPoint.x();
0191   inPar(1) = linPoint.y();
0192   inPar(2) = linPoint.z();
0193
0194   //constraint equations value and number of iterations
0195   double eq;
0196   int nit = 0;
0197   iterations = 0;
0198   csum = 0.0;
0199
0200   GlobalPoint lPoint = linPoint;
0201   RefCountedKinematicVertex rVtx;
0202   ROOT::Math::SMatrix<double, 3 + 7 * nTrk, 3 + 7 * nTrk, ROOT::Math::MatRepSym<double, 3 + 7 * nTrk> > refCCov =
0203       ROOT::Math::SMatrixNoInit();
0204
0205   double chisq = 1e6;
0206   bool convergence = false;
0207
0208   //iterarions over the updator: each time updated parameters
0209   //are taken as new linearization point
0210   do {
0211     eq = 0.;
0212     refCCov = inCov;
0213     std::vector<KinematicState> oldStates = lStates;
0214     GlobalVector mf = field->inInverseGeV(lPoint);
0215     rVtx = updator->update(inPar, refCCov, lStates, lPoint, mf, cs);
0216     if (particles.size() != lStates.size() || rVtx == nullptr) {
0217       LogDebug("KinematicConstrainedVertexFitter") << "updator failure\n";
0218       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0219     }
0220
0221     double newchisq = rVtx->chiSquared();
0222     if (nit > 2 && newchisq > theMaxReducedChiSq * rVtx->degreesOfFreedom() &&
0223         (newchisq - chisq) > (-theMinChiSqImprovement)) {
0224       LogDebug("KinematicConstrainedVertexFitter") << "bad chisq and insufficient improvement, bailing\n";
0225       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0226     }
0227     chisq = newchisq;
0228
0229     const GlobalPoint& newPoint = rVtx->position();
0230
0231     double maxDelta = 0.0;
0232
0233     double deltapos[3];
0234     deltapos[0] = newPoint.x() - lPoint.x();
0235     deltapos[1] = newPoint.y() - lPoint.y();
0236     deltapos[2] = newPoint.z() - lPoint.z();
0237     for (int i = 0; i < 3; ++i) {
0238       double delta = deltapos[i] * deltapos[i] / rVtx->error().matrix()(i, i);
0239       if (delta > maxDelta)
0240         maxDelta = delta;
0241     }
0242
0243     for (std::vector<KinematicState>::const_iterator itold = oldStates.begin(), itnew = lStates.begin();
0244          itnew != lStates.end();
0245          ++itold, ++itnew) {
0246       for (int i = 0; i < 7; ++i) {
0247         double deltapar = itnew->kinematicParameters()(i) - itold->kinematicParameters()(i);
0248         double delta = deltapar * deltapar / itnew->kinematicParametersError().matrix()(i, i);
0249         if (delta > maxDelta)
0250           maxDelta = delta;
0251       }
0252     }
0253
0254     lPoint = newPoint;
0255
0256     nit++;
0257     convergence = maxDelta < theMaxDelta || (nit == theMaxStep && maxDelta < 4.0 * theMaxDelta);
0258
0259   } while (nit < theMaxStep && !convergence);
0260
0261   if (!convergence) {
0262     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0263   }
0264
0265   // std::cout << "new full cov matrix" << std::endl;
0266   // std::cout << refCCov << std::endl;
0267
0268   iterations = nit;
0269   csum = eq;
0270
0271   return tBuilder->buildTree<nTrk>(particles, lStates, rVtx, refCCov);
0272 }
0273
0274 template <int nTrk, int nConstraint>
0275 int KinematicConstrainedVertexFitterT<nTrk, nConstraint>::getNit() const {
0276   return iterations;
0277 }
0278
0279 template <int nTrk, int nConstraint>
0280 float KinematicConstrainedVertexFitterT<nTrk, nConstraint>::getCSum() const {
0281   return csum;
0282 }
0283
0284 #endif