LagrangeParentParticleFitter Class Reference

#include <LagrangeParentParticleFitter.h>

Inheritance diagram for LagrangeParentParticleFitter:

Public Member Functions
LagrangeParentParticleFitter *	clone () const
std::vector< RefCountedKinematicTree >	fit (const std::vector< RefCountedKinematicTree > &trees, KinematicConstraint *cs) const
	LagrangeParentParticleFitter ()
void	setParameters (const edm::ParameterSet &pSet)
	~LagrangeParentParticleFitter ()
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

(

)

Definition at line 6 of file LagrangeParentParticleFitter.cc.

References defaultParameters().

Referenced by clone().

{defaultParameters();}

LagrangeParentParticleFitter::~LagrangeParentParticleFitter ( )

inline

Definition at line 25 of file LagrangeParentParticleFitter.h.

References fwrapper::cs, fit(), and setParameters().

{}

Member Function Documentation

LagrangeParentParticleFitter* LagrangeParentParticleFitter::clone ( ) const

inlinevirtual

Clone method

Implements ParentParticleFitter.

Definition at line 44 of file LagrangeParentParticleFitter.h.

References defaultParameters(), and LagrangeParentParticleFitter().

 {return new LagrangeParentParticleFitter(*this);}

void LagrangeParentParticleFitter::defaultParameters ( )

private

Definition at line 282 of file LagrangeParentParticleFitter.cc.

References theMaxDiff, and theMaxStep.

Referenced by clone(), and LagrangeParentParticleFitter().

{
  theMaxDiff = 0.001;
  theMaxStep = 100;
}

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

virtual

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 121 of file LagrangeParentParticleFitter.cc.

References funct::abs(), KinematicConstraint::clone(), KinematicConstraint::derivative(), KinematicConstraint::deviations(), runTauDisplay::dr, i, j, gen::k, LogDebug, KinematicParametersError::matrix(), KinematicConstraint::numberOfEquations(), InputSort::sort(), theMaxDiff, theMaxStep, KinematicConstraint::value(), and KinematicVertexFactory::vertex().

Referenced by trackingPlots.Iteration::modules(), and ~LagrangeParentParticleFitter().

{
 
 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; 
}     

void LagrangeParentParticleFitter::setParameters

(

const edm::ParameterSet &

pSet

)

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

Definition at line 276 of file LagrangeParentParticleFitter.cc.

References edm::ParameterSet::getParameter(), theMaxDiff, and theMaxStep.

Referenced by ~LagrangeParentParticleFitter().

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

Member Data Documentation

float LagrangeParentParticleFitter::theMaxDiff

private

Definition at line 51 of file LagrangeParentParticleFitter.h.

Referenced by defaultParameters(), fit(), and setParameters().

int LagrangeParentParticleFitter::theMaxStep

private

Definition at line 52 of file LagrangeParentParticleFitter.h.

Referenced by defaultParameters(), fit(), and setParameters().