#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() [1/2]

KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter ( )

Default constructor using LMSLinearizationPointFinder

Definition at line 7 of file KinematicConstrainedVertexFitter.cc.

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

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

◆ KinematicConstrainedVertexFitter() [2/2]

KinematicConstrainedVertexFitter::KinematicConstrainedVertexFitter ( const LinearizationPointFinder & fnd )

Constructor with user-provided LinearizationPointFinder

Definition at line 17 of file KinematicConstrainedVertexFitter.cc.

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

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

◆ ~KinematicConstrainedVertexFitter()

KinematicConstrainedVertexFitter::~KinematicConstrainedVertexFitter ( )

Definition at line 27 of file KinematicConstrainedVertexFitter.cc.

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

References finder, tBuilder, updator, and vCons.

Member Function Documentation

◆ defaultParameters()

void KinematicConstrainedVertexFitter::defaultParameters ( )

private

Definition at line 41 of file KinematicConstrainedVertexFitter.cc.

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

References theMaxDelta, theMaxReducedChiSq, theMaxStep, and theMinChiSqImprovement.

Referenced by KinematicConstrainedVertexFitter().

◆ fit() [1/3]

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 46 of file KinematicConstrainedVertexFitter.h.

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

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

◆ fit() [2/3]

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 53 of file KinematicConstrainedVertexFitter.h.

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

References fwrapper::cs, and fit().

Referenced by trackingPlots.Iteration::modules().

◆ fit() [3/3]

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 48 of file KinematicConstrainedVertexFitter.cc.

                                                                                {
   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 != nullptr) {
     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);
 }

References ConstrainedTreeBuilder::buildTree(), fwrapper::cs, csum, dumpMFGeometry_cfg::delta, edm::eq(), finder, LinearizationPointFinder::getLinearizationPoint(), mps_fire::i, input, KinematicState::isValid(), iterations, dqmiolumiharvest::j, KinematicState::kinematicParameters(), KinematicState::kinematicParametersError(), LogDebug, KinematicParametersError::matrix(), allConversions_cfi::maxDelta, ecalTrigSettings_cff::particles, DiDispStaMuonMonitor_cfi::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().

◆ getCSum()

float KinematicConstrainedVertexFitter::getCSum ( ) const

Definition at line 197 of file KinematicConstrainedVertexFitter.cc.

197 { return csum; }

References csum.

◆ getNit()

int KinematicConstrainedVertexFitter::getNit ( ) const

Definition at line 195 of file KinematicConstrainedVertexFitter.cc.

195 { return iterations; }

References iterations.

Referenced by ConversionVertexFinder::run().

◆ setParameters()

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

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

Definition at line 34 of file KinematicConstrainedVertexFitter.cc.

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