#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 (const LinearizationPointFinder &fnd)
	KinematicConstrainedVertexFitter ()
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 KinematicConstrainedVertexFitter().

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

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 fit().

                                                                {
   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, LaserDQM_cfg::input, KinematicState::isValid(), iterations, j, KinematicState::kinematicParameters(), KinematicState::kinematicParametersError(), LogDebug, KinematicParametersError::matrix(), 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().

{
 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_new()(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);

}

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(), and ConversionVertexFinder::run().

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

float KinematicConstrainedVertexFitter::getCSum ( ) const

Definition at line 209 of file KinematicConstrainedVertexFitter.cc.

References csum.

                                                      {
    return csum;
}

int KinematicConstrainedVertexFitter::getNit ( ) const

Definition at line 205 of file KinematicConstrainedVertexFitter.cc.

References iterations.

Referenced by 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");
}