CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoVertex/KinematicFit/src/MultiTrackPointingKinematicConstraint.cc

Go to the documentation of this file.
00001 #include "RecoVertex/KinematicFit/interface/MultiTrackPointingKinematicConstraint.h"
00002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00003 
00004 AlgebraicVector MultiTrackPointingKinematicConstraint::value(const std::vector<KinematicState> states, const GlobalPoint& point) const{
00005         int num = states.size();
00006         if(num<2) throw VertexException("MultiTrackPointingKinematicConstraint::value <2 states passed");
00007 
00008         //2 equations (for all tracks)
00009         AlgebraicVector  vl(2,0);
00010         double dx = point.x() - refPoint.x();
00011         double dy = point.y() - refPoint.y();
00012         double dz = point.z() - refPoint.z();
00013         double dT = sqrt(pow(dx,2) + pow(dy,2));
00014         double ds = sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2));
00015 
00016         double pxSum=0, pySum=0, pzSum=0;
00017         for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++)
00018         {
00019                 double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
00020 
00021                 pxSum += i->kinematicParameters()(3) - a*(point.y() - i->kinematicParameters()(1));
00022                 pySum += i->kinematicParameters()(4) + a*(point.x() - i->kinematicParameters()(0));
00023                 pzSum += i->kinematicParameters()(5);
00024         }
00025         
00026         double pT = sqrt(pow(pxSum,2) + pow(pySum,2));
00027         double pSum = sqrt(pow(pxSum,2) + pow(pySum,2) + pow(pzSum,2));
00028         
00029         vl(1) = (dT - dx)/dy + (pxSum - pT)/pySum;
00030         vl(2) = (ds - dT)/dz + (pT - pSum)/pzSum;
00031         
00032         return vl;
00033 }
00034 
00035 AlgebraicMatrix MultiTrackPointingKinematicConstraint::parametersDerivative(const std::vector<KinematicState> states, const GlobalPoint& point) const{
00036         int num = states.size();
00037         if(num<2) throw VertexException("MultiTrackPointingKinematicConstraint::parametersDerivative <2 states passed");
00038         
00039         //2 equations (for all tracks)
00040         AlgebraicMatrix  matrix(2,num*7,0);//AlgebraicMatrix starts from 1
00041         
00042         double pxSum=0, pySum=0, pzSum=0;
00043         for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++)
00044         {
00045                 double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
00046                 
00047                 pxSum += i->kinematicParameters()(3) - a*(point.y() - i->kinematicParameters()(1));
00048                 pySum += i->kinematicParameters()(4) + a*(point.x() - i->kinematicParameters()(0));
00049                 pzSum += i->kinematicParameters()(5);
00050         }
00051 
00052         double pT = sqrt(pow(pxSum,2) + pow(pySum,2));
00053         double pSum = sqrt(pow(pxSum,2) + pow(pySum,2) + pow(pzSum,2));
00054 
00055         int col=0;
00056         for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++){
00057                 double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
00058 
00059                 matrix(1,1+col*7) =     a*(1/pT + (-pT + pxSum)/pow(pySum,2));//dH/dx
00060                 matrix(1,2+col*7) =     (a - (a*pxSum)/pT)/pySum;//dH/dy
00061                 //dH/dz=0
00062                 matrix(1,4+col*7) =     (pT - pxSum)/(pT*pySum);//dH/dpx
00063                 matrix(1,5+col*7) =     -(1/pT) + (pT - pxSum)/pow(pySum,2);//dH/dpy            
00064                 //dH/dpz=0
00065                 //dH/dm=0
00066                 matrix(2,1+col*7) =     (a*(-pSum + pT)*pySum)/(pSum*pT*pzSum);//dH/dx
00067                 matrix(2,2+col*7) =     (a*( pSum - pT)*pxSum)/(pSum*pT*pzSum);//dH/dy
00068                 //dH/dz
00069                 matrix(2,4+col*7) =     ((-(1/pSum) + 1/pT)*pxSum)/pzSum;//dH/dpx
00070                 matrix(2,5+col*7) =     ((-(1/pSum) + 1/pT)*pySum)/pzSum;//dH/dpy
00071                 matrix(2,6+col*7) =     -(1/pSum) + (pSum - pT)/pow(pzSum,2);//dH/dpz
00072                 //dH/dm=0               
00073 
00074                 col++;
00075         }
00076         
00077         return matrix;
00078 }
00079 
00080 AlgebraicMatrix MultiTrackPointingKinematicConstraint::positionDerivative(const std::vector<KinematicState> states, const GlobalPoint& point) const{
00081         int num = states.size();
00082         if(num<2) throw VertexException("MultiTrackPointingKinematicConstraint::positionDerivative <2 states passed");
00083         
00084         //2 equations (for all tracks)
00085         AlgebraicMatrix  matrix(2,3,0);
00086         double dx = point.x() - refPoint.x();
00087         double dy = point.y() - refPoint.y();
00088         double dz = point.z() - refPoint.z();
00089         double dT = sqrt(pow(dx,2) + pow(dy,2));
00090         double ds = sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2));
00091         
00092         double pxSum=0, pySum=0, pzSum=0, aSum = 0;
00093         for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++){
00094                 double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
00095                 aSum += a;
00096                 
00097                 pxSum += i->kinematicParameters()(3) - a*(point.y() - i->kinematicParameters()(1));
00098                 pySum += i->kinematicParameters()(4) + a*(point.x() - i->kinematicParameters()(0));
00099                 pzSum += i->kinematicParameters()(5);
00100         }
00101         double pT = sqrt(pow(pxSum,2) + pow(pySum,2));
00102         double pSum = sqrt(pow(pxSum,2) + pow(pySum,2) + pow(pzSum,2));
00103 
00104         matrix(1,1) = (-1 + dx/dT)/dy - aSum/pT + (aSum*(pT - pxSum))/pow(pySum,2);//dH/dxv
00105         matrix(1,2) = 1/dT + (-dT + dx)/pow(dy,2) - (aSum*(pT - pxSum))/(pT*pySum);//dH/dyv
00106         //dH/dzv=0
00107         matrix(2,1) = ((1/ds - 1/dT)*dx)/dz + (aSum*(pSum - pT)*pySum)/(pSum*pT*pzSum);//dH/dxv
00108         matrix(2,2) = ((1/ds - 1/dT)*dy)/dz - (aSum*(pSum - pT)*pxSum)/(pSum*pT*pzSum);//dH/dyv
00109         matrix(2,3) = 1/ds + (-ds + dT)/pow(dz,2);//dH/dzv
00110         
00111         return matrix;
00112 }
00113 
00114 int MultiTrackPointingKinematicConstraint::numberOfEquations() const{
00115         return 2;
00116 }