CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoVertex/KinematicFit/src/KinematicConstrainedVertexFitter.cc

Go to the documentation of this file.
00001 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
00002 #include "RecoVertex/KinematicFit/interface/InputSort.h"
00003 #include "RecoVertex/LinearizationPointFinders/interface/DefaultLinearizationPointFinder.h"
00004 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicVertexFactory.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 
00008 
00009 KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter()
00010 {
00011  finder = new DefaultLinearizationPointFinder();
00012  vCons = new VertexKinematicConstraint();
00013  updator = new KinematicConstrainedVertexUpdator();
00014  tBuilder = new ConstrainedTreeBuilder;
00015  defaultParameters();
00016  iterations = -1;
00017  csum = -1000.0;
00018 }
00019 
00020 KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter(const LinearizationPointFinder& fnd)
00021 {
00022  finder = fnd.clone();
00023  vCons = new VertexKinematicConstraint();
00024  updator = new KinematicConstrainedVertexUpdator();
00025  tBuilder = new ConstrainedTreeBuilder;
00026  defaultParameters();
00027  iterations = -1;
00028  csum = -1000.0;
00029 }
00030 
00031 KinematicConstrainedVertexFitter::~KinematicConstrainedVertexFitter()
00032 {
00033  delete finder;
00034  delete vCons;
00035  delete updator;
00036  delete tBuilder;
00037 }
00038 
00039 void KinematicConstrainedVertexFitter::setParameters(const edm::ParameterSet& pSet)
00040 {
00041   theMaxDelta = pSet.getParameter<double>("maxDelta");
00042   theMaxStep = pSet.getParameter<int>("maxNbrOfIterations");
00043   theMaxReducedChiSq = pSet.getParameter<double>("maxReducedChiSq");
00044   theMinChiSqImprovement = pSet.getParameter<double>("minChiSqImprovement");
00045 }
00046 
00047 void KinematicConstrainedVertexFitter::defaultParameters()
00048 {
00049   theMaxDelta = 0.01;
00050   theMaxStep = 1000;
00051   theMaxReducedChiSq = 225.;
00052   theMinChiSqImprovement = 50.;
00053 }
00054 
00055 RefCountedKinematicTree KinematicConstrainedVertexFitter::fit(std::vector<RefCountedKinematicParticle> part,
00056                                                              MultiTrackKinematicConstraint * cs,
00057                                                              GlobalPoint * pt)
00058 {
00059  if(part.size()<2) throw VertexException("KinematicConstrainedVertexFitter::input states are less than 2");
00060 
00061 //sorting out the input particles
00062  InputSort iSort;
00063  std::pair<std::vector<RefCountedKinematicParticle>, std::vector<FreeTrajectoryState> > input = iSort.sort(part);
00064  const std::vector<RefCountedKinematicParticle> & particles  = input.first;
00065  const std::vector<FreeTrajectoryState> & fStates = input.second;
00066 
00067 // linearization point
00068 // (only compute it using the linearization point finder if no point was passed to the fit function):
00069  GlobalPoint linPoint;
00070  if (pt!=0) {
00071    linPoint  = *pt;
00072  }
00073  else {
00074    linPoint = finder->getLinearizationPoint(fStates);
00075  }
00076 
00077 //initial parameters:
00078  int vSize = particles.size();
00079  AlgebraicVector inPar(3 + 7*vSize,0);
00080 
00081 //final parameters
00082  AlgebraicVector finPar(3 + 7*vSize,0);
00083 
00084 //initial covariance
00085  AlgebraicMatrix inCov(3 + 7*vSize,3 + 7*vSize,0);
00086 
00087 //making initial vector of parameters and initial particle-related covariance
00088  int nSt = 0;
00089  std::vector<KinematicState> inStates;
00090  for(std::vector<RefCountedKinematicParticle>::const_iterator i = particles.begin(); i!=particles.end(); i++)
00091  {
00092   KinematicState state = (*i)->stateAtPoint(linPoint);
00093   if (!state.isValid()) {
00094       LogDebug("KinematicConstrainedVertexFitter")
00095        << "State is invalid at point: "<<linPoint<<std::endl;
00096       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
00097   }
00098   AlgebraicVector prPar = asHepVector<7>(state.kinematicParameters().vector());
00099   for(int j = 1; j<8; j++){inPar(3 + 7*nSt + j) = prPar(j);}
00100   AlgebraicSymMatrix l_cov  = asHepMatrix<7>(state.kinematicParametersError().matrix());
00101   inCov.sub(4 + 7*nSt,4 + 7*nSt ,l_cov);
00102   inStates.push_back(state);
00103   ++nSt;
00104  }
00105 
00106 //initial vertex error matrix components (huge error method)
00107 //and vertex related initial vector components
00108  double in_er = 100.;
00109  inCov(1,1) = in_er;
00110  inCov(2,2) = in_er;
00111  inCov(3,3) = in_er;
00112 
00113  inPar(1) = linPoint.x();
00114  inPar(2) = linPoint.y();
00115  inPar(3) = linPoint.z();
00116 
00117 //constraint equations value and number of iterations
00118  double eq;
00119  int nit = 0;
00120  iterations = 0;
00121  csum = 0.0;
00122 
00123  std::vector<KinematicState> lStates = inStates;
00124  GlobalPoint lPoint  = linPoint;
00125  RefCountedKinematicVertex rVtx;
00126  AlgebraicMatrix refCCov;
00127 
00128  double chisq = 1e6;
00129  bool convergence = false;
00130 //iterarions over the updator: each time updated parameters
00131 //are taken as new linearization point
00132  do{
00133   eq = 0.;
00134   std::pair< std::pair< std::vector<KinematicState>, AlgebraicMatrix >,RefCountedKinematicVertex> lRes =
00135                                       updator->update(inPar,inCov,lStates,lPoint,cs);
00136  
00137   const std::vector<KinematicState> &newStates = lRes.first.first;
00138 
00139   if (particles.size() != newStates.size()) {
00140     LogDebug("KinematicConstrainedVertexFitter")
00141         << "updator failure\n";
00142     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
00143   }
00144 
00145                                       
00146   rVtx = lRes.second;                                      
00147                     
00148   double newchisq = rVtx->chiSquared();
00149   if ( nit>2 && newchisq > theMaxReducedChiSq*rVtx->degreesOfFreedom() && (newchisq-chisq) > (-theMinChiSqImprovement) ) {
00150     LogDebug("KinematicConstrainedVertexFitter")
00151            << "bad chisq and insufficient improvement, bailing\n";
00152     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
00153   }
00154   chisq = newchisq;
00155   
00156 
00157   const GlobalPoint &newPoint = rVtx->position();
00158   
00159   double maxDelta = 0.0;
00160   
00161   double deltapos[3];
00162   deltapos[0] = newPoint.x() - lPoint.x();
00163   deltapos[1] = newPoint.y() - lPoint.y();
00164   deltapos[2] = newPoint.z() - lPoint.z();
00165   for (int i=0; i<3; ++i) {
00166     double delta = deltapos[i]*deltapos[i]/rVtx->error().matrix_new()(i,i);
00167     if (delta>maxDelta) maxDelta = delta;
00168   }
00169   
00170   for (std::vector<KinematicState>::const_iterator itold = lStates.begin(), itnew = newStates.begin();
00171        itnew!=newStates.end(); ++itold,++itnew) {
00172     for (int i=0; i<7; ++i) {
00173       double deltapar = itnew->kinematicParameters()(i) - itold->kinematicParameters()(i);
00174       double delta = deltapar*deltapar/itnew->kinematicParametersError().matrix()(i,i);
00175       if (delta>maxDelta) maxDelta = delta;
00176     }
00177   }
00178   
00179   lStates = newStates;
00180   lPoint = newPoint;
00181 
00182   refCCov = lRes.first.second;
00183   nit++;
00184   convergence = maxDelta<theMaxDelta || (nit==theMaxStep && maxDelta<4.0*theMaxDelta);
00185 
00186  }while(nit<theMaxStep && !convergence);
00187 
00188  if (!convergence) {
00189    return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
00190  } 
00191 
00192   // std::cout << "old full cov matrix" << std::endl;
00193   // std::cout << refCCov << std::endl;
00194 
00195 
00196 // cout<<"number of relinearizations "<<nit<<endl;
00197 // cout<<"value obtained: "<<eq<<endl;
00198   iterations = nit;
00199   csum = eq;
00200 
00201   return  tBuilder->buildTree(particles, lStates, rVtx, refCCov);
00202 
00203 }
00204 
00205 int KinematicConstrainedVertexFitter::getNit() const {
00206     return iterations;
00207 }
00208 
00209 float KinematicConstrainedVertexFitter::getCSum() const {
00210     return csum;
00211 }