Reference Manual

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