KinematicFit/src/KinematicParticleVertexFitter.cc

0001 #include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
0002 // #include "Vertex/LinearizationPointFinders/interface/LMSLinearizationPointFinder.h"
0003 #include "RecoVertex/KinematicFit/interface/FinalTreeBuilder.h"
0004 #include "RecoVertex/VertexTools/interface/SequentialVertexSmoother.h"
0005 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
0006 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackUpdator.h"
0007 #include "RecoVertex/KalmanVertexFit/interface/KalmanSmoothedVertexChi2Estimator.h"
0008 #include "RecoVertex/KalmanVertexFit/interface/KalmanTrackToTrackCovCalculator.h"
0009 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0010 #include "RecoVertex/LinearizationPointFinders/interface/DefaultLinearizationPointFinder.h"
0011 #include "RecoVertex/VertexTools/interface/SequentialVertexFitter.h"
0012 #include "DataFormats/CLHEP/interface/Migration.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 KinematicParticleVertexFitter::KinematicParticleVertexFitter() {
0016   edm::ParameterSet pSet = defaultParameters();
0017   setup(pSet);
0018 }
0019
0020 KinematicParticleVertexFitter::KinematicParticleVertexFitter(const edm::ParameterSet &pSet) { setup(pSet); }
0021
0022 void KinematicParticleVertexFitter::setup(const edm::ParameterSet &pSet) {
0023   pointFinder = new DefaultLinearizationPointFinder();
0024   vFactory = new VertexTrackFactory<6>();
0025
0026   KalmanVertexTrackUpdator<6> vtu;
0027   KalmanSmoothedVertexChi2Estimator<6> vse;
0028   KalmanTrackToTrackCovCalculator<6> covCalc;
0029   SequentialVertexSmoother<6> smoother(vtu, vse, covCalc);
0030   fitter = new SequentialVertexFitter<6>(
0031       pSet, *pointFinder, KalmanVertexUpdator<6>(), smoother, ParticleKinematicLinearizedTrackStateFactory());
0032 }
0033
0034 KinematicParticleVertexFitter::~KinematicParticleVertexFitter() {
0035   delete vFactory;
0036   delete pointFinder;
0037   delete fitter;
0038 }
0039
0040 edm::ParameterSet KinematicParticleVertexFitter::defaultParameters() const {
0041   edm::ParameterSet pSet;
0042   pSet.addParameter<double>("maxDistance", 0.01);
0043   pSet.addParameter<int>("maxNbrOfIterations", 100);  //10
0044   return pSet;
0045 }
0046
0047 RefCountedKinematicTree KinematicParticleVertexFitter::fit(
0048     const std::vector<RefCountedKinematicParticle> &particles) const {
0049   typedef ReferenceCountingPointer<VertexTrack<6> > RefCountedVertexTrack;
0050   //sorting the input
0051   if (particles.size() < 2)
0052     throw VertexException("KinematicParticleVertexFitter::input states are less than 2");
0053   InputSort iSort;
0054   std::pair<std::vector<RefCountedKinematicParticle>, std::vector<FreeTrajectoryState> > input = iSort.sort(particles);
0055   std::vector<RefCountedKinematicParticle> &newPart = input.first;
0056   std::vector<FreeTrajectoryState> &freeStates = input.second;
0057
0058   GlobalPoint linPoint = pointFinder->getLinearizationPoint(freeStates);
0059
0060   // cout<<"Linearization point found"<<endl;
0061
0062   //making initial veretx seed with lin point as position and a fake error
0063   AlgebraicSymMatrix33 we;
0064   we(0, 0) = we(1, 1) = we(2, 2) = 10000.;
0065   GlobalError error(we);
0066   VertexState state(linPoint, error);
0067
0068   //vector of Vertex Tracks to fit
0069   std::vector<RefCountedVertexTrack> ttf;
0070   TrackKinematicStatePropagator propagator_;
0071   for (auto const &i : newPart) {
0072     if (!(i)->currentState().isValid() || !propagator_.willPropagateToTheTransversePCA((i)->currentState(), linPoint)) {
0073       // std::cout << "Here's the bad state." << std::endl;
0074       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());  // return invalid vertex
0075     }
0076     ttf.push_back(vFactory->vertexTrack((i)->particleLinearizedTrackState(linPoint), state, 1.));
0077   }
0078
0079   // //debugging code to check neutrals:
0080   //  for(std::vector<RefCountedVertexTrack>::const_iterator i = ttf.begin(); i!=ttf.end(); i++)
0081   //  {
0082   // //   cout<<"predicted state momentum error"<<(*i)->linearizedTrack()->predictedStateMomentumError()<<endl;
0083   // //  cout<<"Momentum jacobian"<<(*i)->linearizedTrack()->momentumJacobian() <<endl;
0084   //  //  cout<<"predicted state momentum "<<(*i)->linearizedTrack()->predictedStateMomentum()<<endl;
0085   // //   cout<<"constant term"<<(*i)->linearizedTrack()->constantTerm()<<endl;
0086   //
0087   //  }
0088
0089   CachingVertex<6> vtx = fitter->vertex(ttf);
0090   if (!vtx.isValid()) {
0091     LogDebug("RecoVertex/KinematicParticleVertexFitter") << "Fitted position is invalid. Returned Tree is invalid\n";
0092     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());  // return invalid vertex
0093   }
0094   FinalTreeBuilder tBuilder;
0095   return tBuilder.buildTree(vtx, newPart);
0096 }