CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
LagrangeParentParticleFitter Class Reference

#include <LagrangeParentParticleFitter.h>

Inheritance diagram for LagrangeParentParticleFitter:
ParentParticleFitter

Public Member Functions

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

LagrangeParentParticleFitter::~LagrangeParentParticleFitter ( )
inline

Definition at line 25 of file LagrangeParentParticleFitter.h.

25 {}

Member Function Documentation

LagrangeParentParticleFitter* LagrangeParentParticleFitter::clone ( ) const
inlinevirtual

Clone method

Implements ParentParticleFitter.

Definition at line 44 of file LagrangeParentParticleFitter.h.

References LagrangeParentParticleFitter().

void LagrangeParentParticleFitter::defaultParameters ( )
private
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(), i, j, relval_steps::k, LogDebug, KinematicParametersError::matrix(), KinematicConstraint::numberOfEquations(), InputSort::sort(), theMaxDiff, theMaxStep, KinematicConstraint::value(), and KinematicVertexFactory::vertex().

123 {
124 
125  InputSort iSort;
126  std::vector<RefCountedKinematicParticle> prt = iSort.sort(trees);
127  int nStates = prt.size();
128 
129 //making full initial parameters and covariance
130  AlgebraicVector part(7*nStates,0);
131  AlgebraicSymMatrix cov(7*nStates,0);
132 
133  AlgebraicVector chi_in(nStates,0);
134  AlgebraicVector ndf_in(nStates,0);
135  int l_c=0;
136  for(std::vector<RefCountedKinematicParticle>::const_iterator i = prt.begin(); i != prt.end(); i++)
137  {
138  AlgebraicVector7 lp = (*i)->currentState().kinematicParameters().vector();
139  for(int j = 1; j != 8; j++){part(7*l_c + j) = lp(j-1);}
140  AlgebraicSymMatrix lc= asHepMatrix<7>((*i)->currentState().kinematicParametersError().matrix());
141  cov.sub(7*l_c+1,lc);
142  chi_in(l_c+1) = (*i)->chiSquared();
143  ndf_in(l_c+1) = (*i)->degreesOfFreedom();
144  l_c++;
145  }
146 //refitted parameters and covariance matrix:
147 //simple and symmetric
148  AlgebraicVector refPar;
149  AlgebraicSymMatrix refCovS;
150 
151 //constraint values, derivatives and deviations:
152  AlgebraicVector vl;
153  AlgebraicMatrix dr;
154  AlgebraicVector dev;
155  int nstep = 0;
156  double df = 0.;
157  AlgebraicVector exPoint = part;
158 
159 //this piece of code should later be refactored:
160 // The algorithm is the same as above, but the
161 // number of refitted particles is >1. Smart way of
162 // refactoring should be chosen for it.
163  AlgebraicVector chi;
164  AlgebraicVector ndf;
165 // cout<<"Starting the main loop"<<endl;
166  do{
167  df = 0.;
168  chi = chi_in;
169  ndf = ndf_in;
170 // cout<<"Iterational linearization point: "<<exPoint<<endl;
171  vl = cs->value(exPoint).first;
172  dr = cs->derivative(exPoint).first;
173  dev = cs->deviations(nStates);
174 
175 // cout<<"The value : "<<vl<<endl;
176 // cout<<"The derivative: "<<dr<<endl;
177 // cout<<"deviations: "<<dev<<endl;
178 // cout<<"covariance "<<cov<<endl;
179 
180 //residual between expansion and current parameters
181 //0 at the first step
182  AlgebraicVector delta_alpha = part - exPoint;
183 
184 //parameters needed for refit
185 // v_d = (D * V_alpha & * D.T)^(-1)
186  AlgebraicMatrix drt = dr.T();
187  AlgebraicMatrix v_d = dr * cov * drt;
188  int ifail = 0;
189  v_d.invert(ifail);
190  if(ifail != 0) {
191  LogDebug("KinematicConstrainedVertexFitter")
192  << "Fit failed: unable to invert covariance matrix\n";
193  return std::vector<RefCountedKinematicTree>();
194  }
195 
196 //lagrangian multipliers
197 //lambda = V_d * (D * delta_alpha + d)
198  AlgebraicVector lambda = v_d * (dr*delta_alpha + vl);
199 
200 //refitted parameters
201  refPar = part - (cov * drt * lambda);
202 
203 //refitted covariance: simple and SymMatrix
204  refCovS = cov;
205  AlgebraicMatrix sPart = drt * v_d * dr;
206  AlgebraicMatrix covF = cov * sPart * cov;
207 
208 //total covariance symmatrix
209  AlgebraicSymMatrix sCovF(nStates*7,0);
210  for(int i = 1; i< nStates*7 +1; ++i)
211  {
212  for(int j = 1; j< nStates*7 +1; j++)
213  {if(i<=j) sCovF(i,j) = covF(i,j);}
214  }
215 
216  refCovS -= sCovF;
217 
218 // cout<<"Fitter: refitted covariance "<<refCovS<<endl;
219  for(int i = 1; i < nStates+1; i++)
220  {for(int j = 1; j<8; j++){refCovS((i-1)+j,(i-1)+j) += dev(j);}}
221 
222 //chiSquared
223  for(int k =1; k<nStates+1; k++)
224  {
225  chi(k) += (lambda.T() * (dr*delta_alpha + vl))(1);
226  ndf(k) += cs->numberOfEquations();
227  }
228 //new expansionPoint
229  exPoint = refPar;
230  AlgebraicVector vlp = cs->value(exPoint).first;
231  for(int i = 1; i< vl.num_row();i++)
232  {df += std::abs(vlp(i));}
233  nstep++;
234  }while(df>theMaxDiff && nstep<theMaxStep);
235 //here math and iterative part is finished, starting an output production
236 //creating an output KinematicParticle and putting it on its place in the tree
237 
238 //vector of refitted particles and trees
239  std::vector<RefCountedKinematicParticle> refPart;
240  std::vector<RefCountedKinematicTree> refTrees = trees;
241 
242  int j=1;
243  std::vector<RefCountedKinematicTree>::const_iterator tr = refTrees.begin();
244  for(std::vector<RefCountedKinematicParticle>::const_iterator i = prt.begin(); i!= prt.end(); i++)
245  {
246  AlgebraicVector7 lRefPar;
247  for(int k = 1; k<8 ; k++)
248  {lRefPar(k-1) = refPar((j-1)*7+k);}
249  AlgebraicSymMatrix77 lRefCovS = asSMatrix<7>(refCovS.sub((j-1)*7 +1,(j-1)*7+7));
250 
251 //new refitted parameters and covariance
252  KinematicParameters param(lRefPar);
253  KinematicParametersError er(lRefCovS);
254  KinematicState kState(param,er,(*i)->initialState().particleCharge(), (**i).magneticField());
255  RefCountedKinematicParticle refParticle = (*i)->refittedParticle(kState,chi(j),ndf(j),cs->clone());
256 
257 //replacing particle itself
258  (*tr)->findParticle(*i);
259  RefCountedKinematicVertex inVertex = (*tr)->currentDecayVertex();
260  (*tr)->replaceCurrentParticle(refParticle);
261 
262 //replacing the vertex with its refitted version
263  GlobalPoint nvPos(param.position());
264  AlgebraicSymMatrix nvMatrix = asHepMatrix<7>(er.matrix()).sub(1,3);
265  GlobalError nvError(asSMatrix<3>(nvMatrix));
266  VertexState vState(nvPos, nvError, 1.0);
267  KinematicVertexFactory vFactory;
268  RefCountedKinematicVertex nVertex = vFactory.vertex(vState,inVertex,chi(j),ndf(j));
269  (*tr)->replaceCurrentVertex(nVertex);
270  tr++;
271  j++;
272  }
273  return refTrees;
274 }
#define LogDebug(id)
int i
Definition: DBlmapReader.cc:9
virtual KinematicConstraint * clone() const =0
ROOT::Math::SMatrix< double, 7, 7, ROOT::Math::MatRepSym< double, 7 > > AlgebraicSymMatrix77
Definition: Matrices.h:16
std::pair< std::vector< RefCountedKinematicParticle >, std::vector< FreeTrajectoryState > > sort(const std::vector< RefCountedKinematicParticle > &particles) const
Definition: InputSort.cc:6
ROOT::Math::SVector< double, 7 > AlgebraicVector7
Definition: Matrices.h:15
CLHEP::HepMatrix AlgebraicMatrix
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
virtual int numberOfEquations() const =0
CLHEP::HepVector AlgebraicVector
virtual AlgebraicVector deviations(int nStates) const =0
part
Definition: HCALResponse.h:20
virtual std::pair< AlgebraicMatrix, AlgebraicVector > derivative(const AlgebraicVector &exPoint) const =0
CLHEP::HepSymMatrix AlgebraicSymMatrix
static RefCountedKinematicVertex vertex(const VertexState &state, float totalChiSq, float degreesOfFr)
virtual std::pair< AlgebraicVector, AlgebraicVector > value(const AlgebraicVector &exPoint) const =0
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.

277 {
278  theMaxDiff = pSet.getParameter<double>("maxDistance");
279  theMaxStep = pSet.getParameter<int>("maxNbrOfIterations");;
280 }
T getParameter(std::string const &) const

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