#include <KinematicConstrainedVertexFitter.h>

Public Member Functions
RefCountedKinematicTree	fit (const std::vector< RefCountedKinematicParticle > &part)

RefCountedKinematicTree	fit (const std::vector< RefCountedKinematicParticle > &part, MultiTrackKinematicConstraint *cs)

RefCountedKinematicTree	fit (const std::vector< RefCountedKinematicParticle > &part, MultiTrackKinematicConstraint cs, GlobalPoint pt)

float	getCSum () const

int	getNit () const

	KinematicConstrainedVertexFitter ()

	KinematicConstrainedVertexFitter (const LinearizationPointFinder &fnd)

void	setParameters (const edm::ParameterSet &pSet)

	~KinematicConstrainedVertexFitter ()

Private Member Functions
void	defaultParameters ()

Private Attributes
float	csum

LinearizationPointFinder *	finder

int	iterations

ConstrainedTreeBuilder *	tBuilder

float	theMaxDelta

float	theMaxReducedChiSq

int	theMaxStep

float	theMinChiSqImprovement

KinematicConstrainedVertexUpdator *	updator

VertexKinematicConstraint *	vCons

Detailed Description

Class fitting the veretx out of set of tracks via usual LMS with Lagrange multipliers. Additional constraints can be applyed to the tracks during the vertex fit (solves non-factorizabele cases). Since the vertex constraint is included by default, do not add a separate VertexKinematicConstraint! Example: Vertex fit with collinear tracks..

Definition at line 21 of file KinematicConstrainedVertexFitter.h.

Constructor & Destructor Documentation

KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter ( )

Default constructor using LMSLinearizationPointFinder

Definition at line 9 of file KinematicConstrainedVertexFitter.cc.

References csum, defaultParameters(), finder, iterations, tBuilder, updator, and vCons.

 {
  finder = new DefaultLinearizationPointFinder();
  vCons = new VertexKinematicConstraint();
  updator = new KinematicConstrainedVertexUpdator();
  tBuilder = new ConstrainedTreeBuilder;
  defaultParameters();
  iterations = -1;
  csum = -1000.0;
 }

KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter ( const LinearizationPointFinder & fnd )

Constructor with user-provided LinearizationPointFinder

Definition at line 20 of file KinematicConstrainedVertexFitter.cc.

References LinearizationPointFinder::clone(), csum, defaultParameters(), finder, iterations, tBuilder, updator, and vCons.

 {
  finder = fnd.clone();
  vCons = new VertexKinematicConstraint();
  updator = new KinematicConstrainedVertexUpdator();
  tBuilder = new ConstrainedTreeBuilder;
  defaultParameters();
  iterations = -1;
  csum = -1000.0;
 }

KinematicConstrainedVertexFitter::~KinematicConstrainedVertexFitter ( )

Definition at line 31 of file KinematicConstrainedVertexFitter.cc.

References finder, tBuilder, updator, and vCons.

 {
  delete finder;
  delete vCons;
  delete updator;
  delete tBuilder;
 }

Member Function Documentation

void KinematicConstrainedVertexFitter::defaultParameters ( )

private

Definition at line 47 of file KinematicConstrainedVertexFitter.cc.

References theMaxDelta, theMaxReducedChiSq, theMaxStep, and theMinChiSqImprovement.

Referenced by fit(), and KinematicConstrainedVertexFitter().

 {
   theMaxDelta = 0.01;
   theMaxStep = 1000;
   theMaxReducedChiSq = 225.;
   theMinChiSqImprovement = 50.;
 }

RefCountedKinematicTree KinematicConstrainedVertexFitter::fit ( const std::vector< RefCountedKinematicParticle > & part )

inline

Without additional constraint, this will perform a simple vertex fit using LMS with Lagrange multipliers method.

Definition at line 48 of file KinematicConstrainedVertexFitter.h.

Referenced by KineExample::analyze(), fit(), BPHKinematicFit::kinematicTree(), trackingPlots.Iteration::modules(), and ConversionVertexFinder::run().

                                                                                  {
    return fit(part, 0, 0);
  }

RefCountedKinematicTree KinematicConstrainedVertexFitter::fit	(	const std::vector< RefCountedKinematicParticle > &	part,
		MultiTrackKinematicConstraint *	cs
	)

inline

LMS with Lagrange multipliers fit of vertex constraint and user-specified constraint.

Definition at line 55 of file KinematicConstrainedVertexFitter.h.

References fwrapper::cs, defaultParameters(), fit(), getCSum(), getNit(), and EnergyCorrector::pt.

Referenced by trackingPlots.Iteration::modules().

                                                                 {
    return fit(part, cs, 0);
  };

RefCountedKinematicTree KinematicConstrainedVertexFitter::fit	(	const std::vector< RefCountedKinematicParticle > &	part,
		MultiTrackKinematicConstraint *	cs,
		GlobalPoint *	pt
	)

LMS with Lagrange multipliers fit of vertex constraint, user-specified constraint and user-specified starting point.

Definition at line 55 of file KinematicConstrainedVertexFitter.cc.

References ConstrainedTreeBuilder::buildTree(), csum, delta, Vispa.Plugins.EdmBrowser.EdmDataAccessor::eq(), finder, LinearizationPointFinder::getLinearizationPoint(), i, input, KinematicState::isValid(), iterations, j, KinematicState::kinematicParameters(), KinematicState::kinematicParametersError(), LogDebug, KinematicParametersError::matrix(), allConversions_cfi::maxDelta, HadronAndPartonSelector_cfi::particles, EnergyCorrector::pt, InputSort::sort(), tBuilder, theMaxDelta, theMaxReducedChiSq, theMaxStep, theMinChiSqImprovement, KinematicConstrainedVertexUpdator::update(), updator, KinematicParameters::vector(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by trackingPlots.Iteration::modules().

 {
  if(part.size()<2) throw VertexException("KinematicConstrainedVertexFitter::input states are less than 2");
 
 //sorting out the input particles
  InputSort iSort;
  std::pair<std::vector<RefCountedKinematicParticle>, std::vector<FreeTrajectoryState> > input = iSort.sort(part);
  const std::vector<RefCountedKinematicParticle> & particles  = input.first;
  const std::vector<FreeTrajectoryState> & fStates = input.second;
 
 // linearization point
 // (only compute it using the linearization point finder if no point was passed to the fit function):
  GlobalPoint linPoint;
  if (pt!=0) {
    linPoint  = *pt;
  }
  else {
    linPoint = finder->getLinearizationPoint(fStates);
  }
 
 //initial parameters:
  int vSize = particles.size();
  AlgebraicVector inPar(3 + 7*vSize,0);
 
 //final parameters
  AlgebraicVector finPar(3 + 7*vSize,0);
 
 //initial covariance
  AlgebraicMatrix inCov(3 + 7*vSize,3 + 7*vSize,0);
 
 //making initial vector of parameters and initial particle-related covariance
  int nSt = 0;
  std::vector<KinematicState> inStates;
  for(std::vector<RefCountedKinematicParticle>::const_iterator i = particles.begin(); i!=particles.end(); i++)
  {
   KinematicState state = (*i)->stateAtPoint(linPoint);
   if (!state.isValid()) {
       LogDebug("KinematicConstrainedVertexFitter")
        << "State is invalid at point: "<<linPoint<<std::endl;
       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
   }
   AlgebraicVector prPar = asHepVector<7>(state.kinematicParameters().vector());
   for(int j = 1; j<8; j++){inPar(3 + 7*nSt + j) = prPar(j);}
   AlgebraicSymMatrix l_cov  = asHepMatrix<7>(state.kinematicParametersError().matrix());
   inCov.sub(4 + 7*nSt,4 + 7*nSt ,l_cov);
   inStates.push_back(state);
   ++nSt;
  }
 
 //initial vertex error matrix components (huge error method)
 //and vertex related initial vector components
  double in_er = 100.;
  inCov(1,1) = in_er;
  inCov(2,2) = in_er;
  inCov(3,3) = in_er;
 
  inPar(1) = linPoint.x();
  inPar(2) = linPoint.y();
  inPar(3) = linPoint.z();
 
 //constraint equations value and number of iterations
  double eq;
  int nit = 0;
  iterations = 0;
  csum = 0.0;
 
  std::vector<KinematicState> lStates = inStates;
  GlobalPoint lPoint  = linPoint;
  RefCountedKinematicVertex rVtx;
  AlgebraicMatrix refCCov;
 
  double chisq = 1e6;
  bool convergence = false;
 //iterarions over the updator: each time updated parameters
 //are taken as new linearization point
  do{
   eq = 0.;
   std::pair< std::pair< std::vector<KinematicState>, AlgebraicMatrix >,RefCountedKinematicVertex> lRes =
                                       updator->update(inPar,inCov,lStates,lPoint,cs);
  
   const std::vector<KinematicState> &newStates = lRes.first.first;
 
   if (particles.size() != newStates.size()) {
     LogDebug("KinematicConstrainedVertexFitter")
         << "updator failure\n";
     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
   }
 
                                       
   rVtx = lRes.second;                                      
                     
   double newchisq = rVtx->chiSquared();
   if ( nit>2 && newchisq > theMaxReducedChiSq*rVtx->degreesOfFreedom() && (newchisq-chisq) > (-theMinChiSqImprovement) ) {
     LogDebug("KinematicConstrainedVertexFitter")
            << "bad chisq and insufficient improvement, bailing\n";
     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
   }
   chisq = newchisq;
   
 
   const GlobalPoint &newPoint = rVtx->position();
   
   double maxDelta = 0.0;
   
   double deltapos[3];
   deltapos[0] = newPoint.x() - lPoint.x();
   deltapos[1] = newPoint.y() - lPoint.y();
   deltapos[2] = newPoint.z() - lPoint.z();
   for (int i=0; i<3; ++i) {
     double delta = deltapos[i]*deltapos[i]/rVtx->error().matrix()(i,i);
     if (delta>maxDelta) maxDelta = delta;
   }
   
   for (std::vector<KinematicState>::const_iterator itold = lStates.begin(), itnew = newStates.begin();
        itnew!=newStates.end(); ++itold,++itnew) {
     for (int i=0; i<7; ++i) {
       double deltapar = itnew->kinematicParameters()(i) - itold->kinematicParameters()(i);
       double delta = deltapar*deltapar/itnew->kinematicParametersError().matrix()(i,i);
       if (delta>maxDelta) maxDelta = delta;
     }
   }
   
   lStates = newStates;
   lPoint = newPoint;
 
   refCCov = lRes.first.second;
   nit++;
   convergence = maxDelta<theMaxDelta || (nit==theMaxStep && maxDelta<4.0*theMaxDelta);
 
  }while(nit<theMaxStep && !convergence);
 
  if (!convergence) {
    return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
  } 
 
   // std::cout << "old full cov matrix" << std::endl;
   // std::cout << refCCov << std::endl;
 
 
 // cout<<"number of relinearizations "<<nit<<endl;
 // cout<<"value obtained: "<<eq<<endl;
   iterations = nit;
   csum = eq;
 
   return  tBuilder->buildTree(particles, lStates, rVtx, refCCov);
 
 }

float KinematicConstrainedVertexFitter::getCSum ( ) const

Definition at line 209 of file KinematicConstrainedVertexFitter.cc.

References csum.

Referenced by fit().

                                                       {
     return csum;
 }

int KinematicConstrainedVertexFitter::getNit ( ) const

Definition at line 205 of file KinematicConstrainedVertexFitter.cc.

References iterations.

Referenced by fit(), and ConversionVertexFinder::run().

                                                    {
     return iterations;
 }

void KinematicConstrainedVertexFitter::setParameters ( const edm::ParameterSet & pSet )

Configuration through PSet: number of iterations(maxDistance) and stopping condition (maxNbrOfIterations)

Definition at line 39 of file KinematicConstrainedVertexFitter.cc.

References edm::ParameterSet::getParameter(), theMaxDelta, theMaxReducedChiSq, theMaxStep, and theMinChiSqImprovement.

 {
   theMaxDelta = pSet.getParameter<double>("maxDelta");
   theMaxStep = pSet.getParameter<int>("maxNbrOfIterations");
   theMaxReducedChiSq = pSet.getParameter<double>("maxReducedChiSq");
   theMinChiSqImprovement = pSet.getParameter<double>("minChiSqImprovement");
 }