CMS 3D CMS Logo

SimplePointingConstraint Class Reference

Topological constraint making a momentum vector to point to the given location in space. More...

#include <RecoVertex/KinematicFit/interface/SimplePointingConstraint.h>

Inheritance diagram for SimplePointingConstraint:

KinematicConstraint

List of all members.

Public Member Functions

virtual SimplePointingConstraintclone () const
 Clone method.
virtual pair< AlgebraicMatrix,
AlgebraicVector
derivative (const vector< RefCountedKinematicParticle > par) const
 Vector of values and matrix of derivatives calculated using current state parameters as expansion point.
virtual pair< AlgebraicMatrix,
AlgebraicVector
derivative (const AlgebraicVector &exPoint) const
virtual AlgebraicVector deviations (int nStates) const
 Returns vector of sigma squared associated to the KinematicParameters of refitted particles Initial deviations are given by user for the constraining parameters (mass, momentum components etc).
virtual int numberOfEquations () const
 Returns number of constraint equations used for fitting.
 SimplePointingConstraint (const GlobalPoint &ref)
virtual pair< AlgebraicVector,
AlgebraicVector
value (const vector< RefCountedKinematicParticle > par) const
 Methods making value and derivative matrix using current state parameters as expansion 7-point.
virtual pair< AlgebraicVector,
AlgebraicVector
value (const AlgebraicVector &exPoint) const
 Vector of values and matrix of derivatives calculated at given expansion 7xNumberOfStates point.

Private Member Functions

pair< AlgebraicMatrix,
AlgebraicVector
makeDerivative (const AlgebraicVector &exPoint) const
pair< AlgebraicVector,
AlgebraicVector
makeValue (const AlgebraicVector &exPoint) const

Private Attributes

GlobalPoint refPoint


Detailed Description

Topological constraint making a momentum vector to point to the given location in space.

Example: if b-meson momentum is reconstructed at b-meson decay position (secondary vertex), making reconstructed momentum pointing the the primary vertex

Multiple track refit is not supported in current version

Kirill Prokofiev, March 2004 MultiState version: July 2004

Definition at line 21 of file SimplePointingConstraint.h.


Constructor & Destructor Documentation

SimplePointingConstraint::SimplePointingConstraint ( const GlobalPoint ref  )  [inline]

Definition at line 25 of file SimplePointingConstraint.h.

Referenced by clone().

00025                                                   :refPoint(ref)
00026   {}


Member Function Documentation

virtual SimplePointingConstraint* SimplePointingConstraint::clone (  )  const [inline, virtual]

Clone method.

Implements KinematicConstraint.

Definition at line 53 of file SimplePointingConstraint.h.

References SimplePointingConstraint().

00054  {return new SimplePointingConstraint(*this);}

pair< AlgebraicMatrix, AlgebraicVector > SimplePointingConstraint::derivative ( const vector< RefCountedKinematicParticle par  )  const [virtual]

Vector of values and matrix of derivatives calculated using current state parameters as expansion point.

Implements KinematicConstraint.

Definition at line 43 of file SimplePointingConstraint.cc.

References makeDerivative().

00044 {
00045  int nStates = par.size();
00046  if(nStates == 0) throw VertexException("PointingKinematicConstraint::Empty vector of particles passed");
00047  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00048  
00049  AlgebraicMatrix dr(2,7,0);
00050  AlgebraicVector lPoint = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
00051 
00052 //2x7 derivative matrix for given state  
00053  AlgebraicMatrix lDeriv = makeDerivative(lPoint).first;
00054  dr.sub(1,1,lDeriv);
00055 // cout<<"Derivative returned: "<<dr<<endl;
00056 // cout<<"For the value: "<<lPoint<<endl;
00057  return pair<AlgebraicMatrix,AlgebraicVector>(dr,lPoint);
00058 }

pair< AlgebraicMatrix, AlgebraicVector > SimplePointingConstraint::derivative ( const AlgebraicVector exPoint  )  const [virtual]

Implements KinematicConstraint.

Definition at line 25 of file SimplePointingConstraint.cc.

References makeDerivative().

00026 {
00027  if(exPoint.num_row() ==0 ) throw VertexException("PointingKinematicConstraint::value requested for zero Linearization point");
00028 
00029 //security check for extended cartesian parametrization 
00030  int inSize = exPoint.num_row(); 
00031  if((inSize%7) !=0) throw VertexException("PointingKinematicConstraint::linearization point has a wrong dimension");
00032  int nStates = inSize/7;
00033  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00034  AlgebraicVector lPar = exPoint;
00035 
00036 //2x7 derivative matrix for given particle
00037  AlgebraicMatrix lDeriv = makeDerivative(lPar).first;
00038  AlgebraicMatrix dr(2,7,0);
00039  dr.sub(1,1,lDeriv);
00040  return pair<AlgebraicMatrix,AlgebraicVector>(dr,lPar);
00041 }

AlgebraicVector SimplePointingConstraint::deviations ( int  nStates  )  const [virtual]

Returns vector of sigma squared associated to the KinematicParameters of refitted particles Initial deviations are given by user for the constraining parameters (mass, momentum components etc).

In case of multiple states exactly the same values are added to every particle parameters

Implements KinematicConstraint.

Definition at line 75 of file SimplePointingConstraint.cc.

00076 {return AlgebraicVector(7*nStates,0);}

pair< AlgebraicMatrix, AlgebraicVector > SimplePointingConstraint::makeDerivative ( const AlgebraicVector exPoint  )  const [private]

Definition at line 120 of file SimplePointingConstraint.cc.

References funct::pow(), refPoint, funct::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by derivative().

00121 { 
00122  AlgebraicMatrix dr(2,7,0);
00123  AlgebraicVector point = exPoint;
00124  double dx = point(1) - refPoint.x();
00125  double dy = point(2) - refPoint.y();
00126  double dz = point(3) - refPoint.z();
00127  double px = point(4);
00128  double py = point(5); 
00129  double pz = point(6);
00130  
00131 
00132 //half angle solution
00133 //d/dx_i
00134  dr(1,1) = (sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2)))) - 
00135             sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2)))))/2.;
00136                                                                   
00137  dr(1,2) = (((-(pow(dx,2)/pow(pow(dx,2) + pow(dy,2),1.5)) + 1/sqrt(pow(dx,2) + pow(dy,2)))*
00138            (1 - px/sqrt(pow(px,2) + pow(py,2))))/
00139            (2.*sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2))))) - 
00140            ((pow(dx,2)/pow(pow(dx,2) + pow(dy,2),1.5) - 1/sqrt(pow(dx,2) + pow(dy,2)))*
00141            (1 + px/sqrt(pow(px,2) + pow(py,2))))/
00142            (2.*sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2))))))/2.;
00143  
00144  
00145  dr(1,3) = 0;
00146  
00147 //d/dp_i  
00148 //debug: x->p index xhange in denominator
00149  dr(1,4) = (-(dx*dy*(1 - px/sqrt(pow(px,2) + pow(py,2))))/
00150            (2.*pow(pow(dx,2) + pow(dy,2),1.5)*
00151            sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2))))) - 
00152            (dx*dy*(1 + px/sqrt(pow(px,2) + pow(py,2))))/
00153            (2.*pow(pow(dx,2) + pow(dy,2),1.5)*
00154            sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2))))))/2.;
00155         
00156                            
00157  dr(1,5) = (((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*px*py)/
00158            (2.*pow(pow(px,2) + pow(py,2),1.5)*
00159            sqrt((1 + dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 - px/sqrt(pow(px,2) + pow(py,2))))) + 
00160            ((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*px*py)/
00161            (2.*pow(pow(px,2) + pow(py,2),1.5)*
00162            sqrt((1 - dx/sqrt(pow(dx,2) + pow(dy,2)))*(1 + px/sqrt(pow(px,2) + pow(py,2))))))/2.;
00163  
00164  
00165  
00166  dr(1,6) = 0;
00167  dr(1,7) = 0; 
00168 
00169 //2nd equation
00170 //d/dx_i
00171 
00172  dr(2,1) =(((-((dx*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)) + 
00173           dx/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00174           (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00175           (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00176           (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00177           (((dx*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5) - 
00178           dx/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00179           (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00180           (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00181           (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00182            
00183  
00184  dr(2,2) = (((-((dy*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)) + 
00185            dy/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00186            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00187            (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00188            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00189            (((dy*sqrt(pow(dx,2) + pow(dy,2)))/pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5) - 
00190            dy/(sqrt(pow(dx,2) + pow(dy,2))*sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2))))*
00191            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00192            (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00193            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00194            
00195         
00196  dr(2,3) = (-(sqrt(pow(dx,2) + pow(dy,2))*dz*(1 - sqrt(pow(px,2) + pow(py,2))/
00197            sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00198            (2.*pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)*
00199            sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00200            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00201            (sqrt(pow(dx,2) + pow(dy,2))*dz*(1 + 
00202            sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))/
00203            (2.*pow(pow(dx,2) + pow(dy,2) + pow(dz,2),1.5)*
00204            sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00205            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00206            
00207 
00208  
00209 //d/dp_i 
00210 //debug: x->p index xhange in denominator
00211 
00212  dr(2,4) = (((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00213            ((px*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5) - 
00214            px/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00215            (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00216            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00217            ((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00218            (-((px*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)) + 
00219            px/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00220            (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00221            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00222            
00223  
00224  dr(2,5) = (((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00225            ((py*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5) - 
00226            py/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00227            (2.*sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00228            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - 
00229            ((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00230            (-((py*sqrt(pow(px,2) + pow(py,2)))/pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)) + 
00231            py/(sqrt(pow(px,2) + pow(py,2))*sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))/
00232            (2.*sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00233            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00234            
00235  
00236  dr(2,6) = (((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00237            sqrt(pow(px,2) + pow(py,2))*pz)/
00238            (2.*pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)*
00239            sqrt((1 + sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00240            (1 - sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + 
00241            ((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00242             sqrt(pow(px,2) + pow(py,2))*pz)/
00243            (2.*pow(pow(px,2) + pow(py,2) + pow(pz,2),1.5)*
00244            sqrt((1 - sqrt(pow(dx,2) + pow(dy,2))/sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2)))*
00245            (1 + sqrt(pow(px,2) + pow(py,2))/sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))))/2.;
00246            
00247  
00248  dr(2,7) = 0;
00249  
00250 // cout<<"derivative matrix "<<dr<<endl;
00251  return pair<AlgebraicMatrix,AlgebraicVector>(dr,point); 
00252 }

pair< AlgebraicVector, AlgebraicVector > SimplePointingConstraint::makeValue ( const AlgebraicVector exPoint  )  const [private]

Definition at line 81 of file SimplePointingConstraint.cc.

References refPoint, funct::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by value().

00082 { 
00083 // cout<<"Make value called"<<endl;
00084  AlgebraicVector vl(2,0);
00085  AlgebraicVector point = exPoint;
00086  double dx = point(1) - refPoint.x();
00087  double dy = point(2) - refPoint.y();
00088  double dz = point(3) - refPoint.z();
00089  double px = point(4);
00090  double py = point(5); 
00091  double pz = point(6);
00092 
00093 
00094 //half angle solution: sin((alpha - betha)/2)
00095  double cos_phi_p = px/sqrt(px*px + py*py);
00096  double cos_phi_x = dx/sqrt(dx*dx + dy*dy);
00097 // cout<<"mom cos phi"<<cos_phi_p<<endl;
00098 // cout<<"x cos phi"<<cos_phi_x<<endl;
00099  
00100  
00101  double cos_theta_p = sqrt(px*px + py*py)/sqrt(px*px + py*py + pz*pz); 
00102  double cos_theta_x = sqrt(dx*dx + dy*dy)/sqrt(dx*dx + dy*dy + dz*dz);
00103  
00104  float feq = sqrt((1-cos_phi_p)*(1+cos_phi_x)) - sqrt((1+cos_phi_p)*(1-cos_phi_x));
00105  float seq = sqrt((1-cos_theta_p)*(1+cos_theta_x)) - sqrt((1+cos_theta_p)*(1-cos_theta_x));
00106  
00107 // cout<<"First factor: "<<feq/2<<endl;
00108 // cout<<"Second factor: "<<seq/2<<endl;
00109  
00110  vl(1) = feq/2;
00111  vl(2) = seq/2;
00112 
00113 // cout<<"Value "<<vl<<endl;
00114 //half angle corrected
00115 // vl(1) = (sin_x/(1+cos_x)) - (sin_p/(1+cos_p));
00116 // vl(2) = (sin_xt/(1+cos_xt)) - (sin_pt/(1+cos_pt));
00117  return pair<AlgebraicVector,AlgebraicVector>(vl,point);
00118 }

int SimplePointingConstraint::numberOfEquations (  )  const [virtual]

Returns number of constraint equations used for fitting.

Method is relevant for proper NDF calculations.

Implements KinematicConstraint.

Definition at line 78 of file SimplePointingConstraint.cc.

00079 {return 2;}

pair< AlgebraicVector, AlgebraicVector > SimplePointingConstraint::value ( const vector< RefCountedKinematicParticle par  )  const [virtual]

Methods making value and derivative matrix using current state parameters as expansion 7-point.

Constraint can be made equaly for single and multiple states

Implements KinematicConstraint.

Definition at line 60 of file SimplePointingConstraint.cc.

References makeValue().

00061 { 
00062  int nStates = par.size();
00063  if(nStates == 0) throw VertexException("PointingKinematicConstraint::Empty vector of particles passed");
00064  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00065  AlgebraicVector vl(2,0);
00066  AlgebraicVector lPoint = asHepVector<7>(par.front()->currentState().kinematicParameters().vector());
00067  vl(1) = makeValue(lPoint).first(1);
00068  vl(2) = makeValue(lPoint).first(2);
00069 // cout<<"Value returned: "<<vl<<endl;
00070 // cout<<"For the point: "<<lPoint<<endl;
00071  
00072  return pair<AlgebraicVector,AlgebraicVector>(vl,lPoint);
00073 }

pair< AlgebraicVector, AlgebraicVector > SimplePointingConstraint::value ( const AlgebraicVector exPoint  )  const [virtual]

Vector of values and matrix of derivatives calculated at given expansion 7xNumberOfStates point.

Implements KinematicConstraint.

Definition at line 5 of file SimplePointingConstraint.cc.

References makeValue().

00006 {
00007  if(exPoint.num_row() ==0 ) throw VertexException("PointingKinematicConstraint::value requested for zero Linearization point");
00008 
00009 //security check for extended cartesian parametrization 
00010  int inSize = exPoint.num_row(); 
00011  if((inSize%7) !=0) throw VertexException("PointingKinematicConstraint::linearization point has a wrong dimension");
00012  int nStates = inSize/7;
00013  if(nStates != 1) throw VertexException("PointingKinematicConstraint::Current version does not support the multistate refit");
00014  
00015  AlgebraicVector lPar = exPoint;
00016  AlgebraicVector vl(2,0);
00017  
00018 //vector of values 1x2  for given particle
00019  AlgebraicVector lValue = makeValue(lPar).first;
00020  vl(1) =lValue(1);
00021  vl(2) =lValue(2);
00022  return pair<AlgebraicVector,AlgebraicVector>(vl,lPar); 
00023 }


Member Data Documentation

GlobalPoint SimplePointingConstraint::refPoint [private]

Definition at line 61 of file SimplePointingConstraint.h.

Referenced by makeDerivative(), and makeValue().


The documentation for this class was generated from the following files:
Generated on Tue Jun 9 18:31:43 2009 for CMSSW by  doxygen 1.5.4