#include <KinematicConstrainedVertexFitterT.h>

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

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

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

float	getCSum () const

int	getNit () const

	KinematicConstrainedVertexFitterT (const MagneticField *ifield)

	KinematicConstrainedVertexFitterT (const MagneticField *ifield, const LinearizationPointFinder &fnd)

void	setParameters (const edm::ParameterSet &pSet)

	~KinematicConstrainedVertexFitterT ()

Private Member Functions
void	defaultParameters ()

Private Attributes
float	csum

const MagneticField *	field

LinearizationPointFinder *	finder

int	iterations

ConstrainedTreeBuilderT *	tBuilder

float	theMaxDelta

float	theMaxReducedChiSq

int	theMaxStep

float	theMinChiSqImprovement

KinematicConstrainedVertexUpdatorT< nTrk, nConstraint > *	updator

VertexKinematicConstraintT *	vCons

Detailed Description

template<int nTrk, int nConstraint>
class KinematicConstrainedVertexFitterT< nTrk, nConstraint >

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 24 of file KinematicConstrainedVertexFitterT.h.

Constructor & Destructor Documentation

template<int nTrk, int nConstraint>

KinematicConstrainedVertexFitterT< nTrk, nConstraint >::KinematicConstrainedVertexFitterT ( const MagneticField * ifield )

explicit

Default constructor using LMSLinearizationPointFinder

Definition at line 104 of file KinematicConstrainedVertexFitterT.h.

References KinematicConstrainedVertexFitterT< nTrk, nConstraint >::csum, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::defaultParameters(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::finder, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::iterations, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::tBuilder, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::updator, and KinematicConstrainedVertexFitterT< nTrk, nConstraint >::vCons.

                                                                                                                     : 
   field(ifield)
 {
   finder = new DefaultLinearizationPointFinder();
   vCons = new VertexKinematicConstraintT();
   updator = new KinematicConstrainedVertexUpdatorT<nTrk,nConstraint>();
   tBuilder = new ConstrainedTreeBuilderT;
   defaultParameters();
   iterations = -1;
   csum = -1000.0;
 }

template<int nTrk, int nConstraint>

KinematicConstrainedVertexFitterT< nTrk, nConstraint >::KinematicConstrainedVertexFitterT	(	const MagneticField *	ifield,
		const LinearizationPointFinder &	fnd
	)

Constructor with user-provided LinearizationPointFinder

Definition at line 117 of file KinematicConstrainedVertexFitterT.h.

References LinearizationPointFinder::clone(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::csum, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::defaultParameters(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::finder, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::iterations, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::tBuilder, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::updator, and KinematicConstrainedVertexFitterT< nTrk, nConstraint >::vCons.

                                                                                                                                                          : 
   field(ifield)
 {
   finder = fnd.clone();
   vCons = new VertexKinematicConstraintT();
   updator = new KinematicConstrainedVertexUpdatorT<nTrk,nConstraint>();
   tBuilder = new ConstrainedTreeBuilderT;
   defaultParameters();
   iterations = -1;
   csum = -1000.0;
 }

template<int nTrk, int nConstraint>

KinematicConstrainedVertexFitterT< nTrk, nConstraint >::~KinematicConstrainedVertexFitterT ( )

Definition at line 130 of file KinematicConstrainedVertexFitterT.h.

References KinematicConstrainedVertexFitterT< nTrk, nConstraint >::finder, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::tBuilder, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::updator, and KinematicConstrainedVertexFitterT< nTrk, nConstraint >::vCons.

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

Member Function Documentation

template<int nTrk, int nConstraint>

void KinematicConstrainedVertexFitterT< nTrk, nConstraint >::defaultParameters ( )

private

Definition at line 148 of file KinematicConstrainedVertexFitterT.h.

References KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMaxDelta, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMaxReducedChiSq, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMaxStep, and KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMinChiSqImprovement.

Referenced by KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit(), and KinematicConstrainedVertexFitterT< nTrk, nConstraint >::KinematicConstrainedVertexFitterT().

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

template<int nTrk, int nConstraint>

RefCountedKinematicTree KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit ( const std::vector< RefCountedKinematicParticle > & part )

inline

Without additional constraint, this will perform a simple vertex fit using LMS with Lagrange multipliers method (by definition valid only if nConstraint=0)

Definition at line 51 of file KinematicConstrainedVertexFitterT.h.

Referenced by KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit(), trackingPlots.Iteration::modules(), and ConversionVertexFinder::run().

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

template<int nTrk, int nConstraint>

RefCountedKinematicTree KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit	(	const std::vector< RefCountedKinematicParticle > &	part,
		MultiTrackKinematicConstraintT< nTrk, nConstraint > *	cs
	)

inline

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

Definition at line 58 of file KinematicConstrainedVertexFitterT.h.

References fwrapper::cs, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::defaultParameters(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::getCSum(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::getNit(), and EnergyCorrector::pt.

Referenced by trackingPlots.Iteration::modules().

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

template<int nTrk, int nConstraint>

RefCountedKinematicTree KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit	(	const std::vector< RefCountedKinematicParticle > &	part,
		MultiTrackKinematicConstraintT< nTrk, nConstraint > *	cs,
		GlobalPoint *	pt
	)

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

Definition at line 159 of file KinematicConstrainedVertexFitterT.h.

Referenced by trackingPlots.Iteration::modules().

 {
    assert( nConstraint==0 || cs!=0);
    if(part.size()!=nTrk) throw VertexException("KinematicConstrainedVertexFitterT::input states are not nTrk");
   
   //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:
   GlobalPoint linPoint  = (pt!=0) ? *pt :  finder->getLinearizationPoint(fStates);
   
   //initial parameters:
   ROOT::Math::SVector<double,3+7*nTrk> inPar; //3+ 7*ntracks
   ROOT::Math::SVector<double,3+7*nTrk> finPar; //3+ 7*ntracks
   
   ROOT::Math::SMatrix<double, 3+7*nTrk,3+7*nTrk ,ROOT::Math::MatRepSym<double,3+7*nTrk> > inCov;
   
   //making initial vector of parameters and initial particle-related covariance
   int nSt = 0;
   std::vector<KinematicState> lStates(nTrk);
   for(std::vector<RefCountedKinematicParticle>::const_iterator i = particles.begin(); i!=particles.end(); i++)
     {
       lStates[nSt] = (*i)->stateAtPoint(linPoint);
       KinematicState const & state = lStates[nSt];
       if (!state.isValid()) {
     LogDebug("KinematicConstrainedVertexFitter")
       << "State is invalid at point: "<<linPoint<<std::endl;
     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
       }
       inPar.Place_at(state.kinematicParameters().vector(),3+7*nSt);
       inCov.Place_at(state.kinematicParametersError().matrix(),3 + 7*nSt,3 + 7*nSt);
       ++nSt;
     }
   
   //initial vertex error matrix components (huge error method)
   //and vertex related initial vector components
   double in_er = 100.;
   inCov(0,0) = in_er;
   inCov(1,1) = in_er;
   inCov(2,2) = in_er;
   
   inPar(0) = linPoint.x();
   inPar(1) = linPoint.y();
   inPar(2) = linPoint.z();
   
   //constraint equations value and number of iterations
   double eq;
   int nit = 0;
   iterations = 0;
   csum = 0.0;
   
   GlobalPoint lPoint  = linPoint;
   RefCountedKinematicVertex rVtx;
   ROOT::Math::SMatrix<double, 3+7*nTrk,3+7*nTrk ,ROOT::Math::MatRepSym<double,3+7*nTrk> > refCCov  = ROOT::Math::SMatrixNoInit();
   
   double chisq = 1e6;
   bool convergence = false;
   
   //iterarions over the updator: each time updated parameters
   //are taken as new linearization point
   do{
     eq = 0.;
     refCCov = inCov;
     std::vector<KinematicState> oldStates = lStates;
     GlobalVector mf = field->inInverseGeV(lPoint);
     rVtx = updator->update(inPar,refCCov,lStates,lPoint,mf,cs);
     if (particles.size() != lStates.size() || rVtx == 0) {
       LogDebug("KinematicConstrainedVertexFitter")
     << "updator failure\n";
       return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
     }
     
     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 = oldStates.begin(), itnew = lStates.begin();
      itnew!=lStates.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;
       }
     }
     
     lPoint = newPoint;
     
     nit++;
     convergence = maxDelta<theMaxDelta || (nit==theMaxStep && maxDelta<4.0*theMaxDelta);
     
   }while(nit<theMaxStep && !convergence);
 
   if (!convergence) {
     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
   } 
   
   // std::cout << "new full cov matrix" << std::endl;
   // std::cout << refCCov << std::endl;  
   
   iterations = nit;
   csum = eq;
   
   return  tBuilder->buildTree<nTrk>(particles, lStates, rVtx, refCCov);
   
 }

template<int nTrk, int nConstraint>

float KinematicConstrainedVertexFitterT< nTrk, nConstraint >::getCSum ( ) const

Definition at line 295 of file KinematicConstrainedVertexFitterT.h.

References KinematicConstrainedVertexFitterT< nTrk, nConstraint >::csum.

Referenced by KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit().

                                                                            {
   return csum;
 }

template<int nTrk, int nConstraint>

int KinematicConstrainedVertexFitterT< nTrk, nConstraint >::getNit ( ) const

Definition at line 290 of file KinematicConstrainedVertexFitterT.h.

References KinematicConstrainedVertexFitterT< nTrk, nConstraint >::iterations.

Referenced by KinematicConstrainedVertexFitterT< nTrk, nConstraint >::fit(), and ConversionVertexFinder::run().

                                                                         {
   return iterations;
 }

template<int nTrk, int nConstraint>

void KinematicConstrainedVertexFitterT< nTrk, nConstraint >::setParameters ( const edm::ParameterSet & pSet )

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

Definition at line 139 of file KinematicConstrainedVertexFitterT.h.

References edm::ParameterSet::getParameter(), KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMaxDelta, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMaxReducedChiSq, KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMaxStep, and KinematicConstrainedVertexFitterT< nTrk, nConstraint >::theMinChiSqImprovement.

Referenced by ConversionVertexFinder::run().

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