#include <LagrangeParentParticleFitter.h>

Inheritance diagram for LagrangeParentParticleFitter:

Public Member Functions
LagrangeParentParticleFitter *	clone () const override

std::vector< RefCountedKinematicTree >	fit (const std::vector< RefCountedKinematicTree > &trees, KinematicConstraint *cs) const override

	LagrangeParentParticleFitter ()

void	setParameters (const edm::ParameterSet &pSet)

	~LagrangeParentParticleFitter () override

Public Member Functions inherited from ParentParticleFitter
	ParentParticleFitter ()

virtual	~ParentParticleFitter ()

Private Member Functions
void	defaultParameters ()

Private Attributes
float	theMaxDiff

int	theMaxStep

Detailed Description

KinematicParticle refit via LMS minimization with Lagrange multipliers method. Tunable for number of iterations and stopping condition. Every new iteration takes the result of previous iteration as a linearization point. This is only needed in case of nonlinear constraints. In the linear case result will be reached after the first step. Stopping condition at the moment: sum of abs. constraint equations values at given point

Definition at line 18 of file LagrangeParentParticleFitter.h.

Constructor & Destructor Documentation

◆ LagrangeParentParticleFitter()

LagrangeParentParticleFitter::LagrangeParentParticleFitter ( )

Definition at line 6 of file LagrangeParentParticleFitter.cc.

6 { defaultParameters(); }

References defaultParameters().

Referenced by clone().

◆ ~LagrangeParentParticleFitter()

LagrangeParentParticleFitter::~LagrangeParentParticleFitter ( )

inlineoverride

Definition at line 22 of file LagrangeParentParticleFitter.h.

22 {}

Member Function Documentation

◆ clone()

LagrangeParentParticleFitter* LagrangeParentParticleFitter::clone ( ) const

inlineoverridevirtual

Clone method

Implements ParentParticleFitter.

Definition at line 41 of file LagrangeParentParticleFitter.h.

41 { return new LagrangeParentParticleFitter(*this); }

References LagrangeParentParticleFitter().

◆ defaultParameters()

void LagrangeParentParticleFitter::defaultParameters ( )

private

Definition at line 283 of file LagrangeParentParticleFitter.cc.

                                                      {
   theMaxDiff = 0.001;
   theMaxStep = 100;
 }

References theMaxDiff, and theMaxStep.

Referenced by LagrangeParentParticleFitter().

◆ fit()

std::vector< RefCountedKinematicTree > LagrangeParentParticleFitter::fit	(	const std::vector< RefCountedKinematicTree > &	trees,
		KinematicConstraint *	cs
	)		const

overridevirtual

Refit method taking KinematicConstraint and vector of trees as imput. Only top particles of corresponding trees are refitted, vertex is not created. Number of trees should be one (single track refit) or greater (multiple track refit). Some constraints may not work with single tracks (back to back for ex.)

Implements ParentParticleFitter.

Definition at line 120 of file LagrangeParentParticleFitter.cc.

                                                                                     {
   InputSort iSort;
   std::vector<RefCountedKinematicParticle> prt = iSort.sort(trees);
   int nStates = prt.size();
  
   //making full initial parameters and covariance
   AlgebraicVector part(7 * nStates, 0);
   AlgebraicSymMatrix cov(7 * nStates, 0);
  
   AlgebraicVector chi_in(nStates, 0);
   AlgebraicVector ndf_in(nStates, 0);
   int l_c = 0;
   for (std::vector<RefCountedKinematicParticle>::const_iterator i = prt.begin(); i != prt.end(); i++) {
     AlgebraicVector7 lp = (*i)->currentState().kinematicParameters().vector();
     for (int j = 1; j != 8; j++) {
       part(7 * l_c + j) = lp(j - 1);
     }
     AlgebraicSymMatrix lc = asHepMatrix<7>((*i)->currentState().kinematicParametersError().matrix());
     cov.sub(7 * l_c + 1, lc);
     chi_in(l_c + 1) = (*i)->chiSquared();
     ndf_in(l_c + 1) = (*i)->degreesOfFreedom();
     l_c++;
   }
   //refitted parameters and covariance matrix:
   //simple and symmetric
   AlgebraicVector refPar;
   AlgebraicSymMatrix refCovS;
  
   //constraint values, derivatives and deviations:
   AlgebraicVector vl;
   AlgebraicMatrix dr;
   AlgebraicVector dev;
   int nstep = 0;
   double df = 0.;
   AlgebraicVector exPoint = part;
  
   //this piece of code should later be refactored:
   // The algorithm is the same as above, but the
   // number of refitted particles is >1. Smart way of
   // refactoring should be chosen for it.
   AlgebraicVector chi;
   AlgebraicVector ndf;
   // cout<<"Starting the main loop"<<endl;
   do {
     df = 0.;
     chi = chi_in;
     ndf = ndf_in;
     //  cout<<"Iterational linearization point: "<<exPoint<<endl;
     vl = cs->value(exPoint).first;
     dr = cs->derivative(exPoint).first;
     dev = cs->deviations(nStates);
  
     //  cout<<"The value : "<<vl<<endl;
     //  cout<<"The derivative: "<<dr<<endl;
     //  cout<<"deviations: "<<dev<<endl;
     //  cout<<"covariance "<<cov<<endl;
  
     //residual between expansion and current parameters
     //0 at the first step
     AlgebraicVector delta_alpha = part - exPoint;
  
     //parameters needed for refit
     // v_d = (D * V_alpha & * D.T)^(-1)
     AlgebraicMatrix drt = dr.T();
     AlgebraicMatrix v_d = dr * cov * drt;
     int ifail = 0;
     v_d.invert(ifail);
     if (ifail != 0) {
       LogDebug("KinematicConstrainedVertexFitter") << "Fit failed: unable to invert covariance matrix\n";
       return std::vector<RefCountedKinematicTree>();
     }
  
     //lagrangian multipliers
     //lambda = V_d * (D * delta_alpha + d)
     AlgebraicVector lambda = v_d * (dr * delta_alpha + vl);
  
     //refitted parameters
     refPar = part - (cov * drt * lambda);
  
     //refitted covariance: simple and SymMatrix
     refCovS = cov;
     AlgebraicMatrix sPart = drt * v_d * dr;
     AlgebraicMatrix covF = cov * sPart * cov;
  
     //total covariance symmatrix
     AlgebraicSymMatrix sCovF(nStates * 7, 0);
     for (int i = 1; i < nStates * 7 + 1; ++i) {
       for (int j = 1; j < nStates * 7 + 1; j++) {
         if (i <= j)
           sCovF(i, j) = covF(i, j);
       }
     }
  
     refCovS -= sCovF;
  
     //  cout<<"Fitter: refitted covariance "<<refCovS<<endl;
     for (int i = 1; i < nStates + 1; i++) {
       for (int j = 1; j < 8; j++) {
         refCovS((i - 1) + j, (i - 1) + j) += dev(j);
       }
     }
  
     //chiSquared
     for (int k = 1; k < nStates + 1; k++) {
       chi(k) += (lambda.T() * (dr * delta_alpha + vl))(1);
       ndf(k) += cs->numberOfEquations();
     }
     //new expansionPoint
     exPoint = refPar;
     AlgebraicVector vlp = cs->value(exPoint).first;
     for (int i = 1; i < vl.num_row(); i++) {
       df += std::abs(vlp(i));
     }
     nstep++;
   } while (df > theMaxDiff && nstep < theMaxStep);
   //here math and iterative part is finished, starting an output production
   //creating an output KinematicParticle and putting it on its place in the tree
  
   //vector of refitted particles and trees
   std::vector<RefCountedKinematicParticle> refPart;
   std::vector<RefCountedKinematicTree> refTrees = trees;
  
   int j = 1;
   std::vector<RefCountedKinematicTree>::const_iterator tr = refTrees.begin();
   for (std::vector<RefCountedKinematicParticle>::const_iterator i = prt.begin(); i != prt.end(); i++) {
     AlgebraicVector7 lRefPar;
     for (int k = 1; k < 8; k++) {
       lRefPar(k - 1) = refPar((j - 1) * 7 + k);
     }
     AlgebraicSymMatrix77 lRefCovS = asSMatrix<7>(refCovS.sub((j - 1) * 7 + 1, (j - 1) * 7 + 7));
  
     //new refitted parameters and covariance
     KinematicParameters param(lRefPar);
     KinematicParametersError er(lRefCovS);
     KinematicState kState(param, er, (*i)->initialState().particleCharge(), (**i).magneticField());
     RefCountedKinematicParticle refParticle = (*i)->refittedParticle(kState, chi(j), ndf(j), cs->clone());
  
     //replacing particle itself
     (*tr)->findParticle(*i);
     RefCountedKinematicVertex inVertex = (*tr)->currentDecayVertex();
     (*tr)->replaceCurrentParticle(refParticle);
  
     //replacing the  vertex with its refitted version
     GlobalPoint nvPos(param.position());
     AlgebraicSymMatrix nvMatrix = asHepMatrix<7>(er.matrix()).sub(1, 3);
     GlobalError nvError(asSMatrix<3>(nvMatrix));
     VertexState vState(nvPos, nvError, 1.0);
     KinematicVertexFactory vFactory;
     RefCountedKinematicVertex nVertex = vFactory.vertex(vState, inVertex, chi(j), ndf(j));
     (*tr)->replaceCurrentVertex(nVertex);
     tr++;
     j++;
   }
   return refTrees;
 }

References funct::abs(), fwrapper::cs, flavorHistoryFilter_cfi::dr, mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, LogDebug, KinematicParametersError::matrix(), L1TEGammaOffline_cfi::nVertex, InputSort::sort(), theMaxDiff, theMaxStep, dqmdumpme::trees, and KinematicVertexFactory::vertex().

Referenced by trackingPlots.Iteration::modules().

◆ setParameters()

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

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

Definition at line 277 of file LagrangeParentParticleFitter.cc.

                                                                             {
   theMaxDiff = pSet.getParameter<double>("maxDistance");
   theMaxStep = pSet.getParameter<int>("maxNbrOfIterations");
   ;
 }