/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoVertex/KinematicFit/src/MultiTrackPointingKinematicConstraint.cc

Go to the documentation of this file.

#include "RecoVertex/KinematicFit/interface/MultiTrackPointingKinematicConstraint.h"
#include "RecoVertex/VertexPrimitives/interface/VertexException.h"

AlgebraicVector MultiTrackPointingKinematicConstraint::value(const std::vector<KinematicState> states, const GlobalPoint& point) const{
        int num = states.size();
        if(num<2) throw VertexException("MultiTrackPointingKinematicConstraint::value <2 states passed");

        //2 equations (for all tracks)
        AlgebraicVector  vl(2,0);
        double dx = point.x() - refPoint.x();
        double dy = point.y() - refPoint.y();
        double dz = point.z() - refPoint.z();
        double dT = sqrt(pow(dx,2) + pow(dy,2));
        double ds = sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2));

        double pxSum=0, pySum=0, pzSum=0;
        for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++)
        {
                double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();

                pxSum += i->kinematicParameters()(3) - a*(point.y() - i->kinematicParameters()(1));
                pySum += i->kinematicParameters()(4) + a*(point.x() - i->kinematicParameters()(0));
                pzSum += i->kinematicParameters()(5);
        }
        
        double pT = sqrt(pow(pxSum,2) + pow(pySum,2));
        double pSum = sqrt(pow(pxSum,2) + pow(pySum,2) + pow(pzSum,2));
        
        vl(1) = (dT - dx)/dy + (pxSum - pT)/pySum;
        vl(2) = (ds - dT)/dz + (pT - pSum)/pzSum;
        
        return vl;
}

AlgebraicMatrix MultiTrackPointingKinematicConstraint::parametersDerivative(const std::vector<KinematicState> states, const GlobalPoint& point) const{
        int num = states.size();
        if(num<2) throw VertexException("MultiTrackPointingKinematicConstraint::parametersDerivative <2 states passed");
        
        //2 equations (for all tracks)
        AlgebraicMatrix  matrix(2,num*7,0);//AlgebraicMatrix starts from 1
        
        double pxSum=0, pySum=0, pzSum=0;
        for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++)
        {
                double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
                
                pxSum += i->kinematicParameters()(3) - a*(point.y() - i->kinematicParameters()(1));
                pySum += i->kinematicParameters()(4) + a*(point.x() - i->kinematicParameters()(0));
                pzSum += i->kinematicParameters()(5);
        }

        double pT = sqrt(pow(pxSum,2) + pow(pySum,2));
        double pSum = sqrt(pow(pxSum,2) + pow(pySum,2) + pow(pzSum,2));

        int col=0;
        for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++){
                double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();

                matrix(1,1+col*7) =     a*(1/pT + (-pT + pxSum)/pow(pySum,2));//dH/dx
                matrix(1,2+col*7) =     (a - (a*pxSum)/pT)/pySum;//dH/dy
                //dH/dz=0
                matrix(1,4+col*7) =     (pT - pxSum)/(pT*pySum);//dH/dpx
                matrix(1,5+col*7) =     -(1/pT) + (pT - pxSum)/pow(pySum,2);//dH/dpy            
                //dH/dpz=0
                //dH/dm=0
                matrix(2,1+col*7) =     (a*(-pSum + pT)*pySum)/(pSum*pT*pzSum);//dH/dx
                matrix(2,2+col*7) =     (a*( pSum - pT)*pxSum)/(pSum*pT*pzSum);//dH/dy
                //dH/dz
                matrix(2,4+col*7) =     ((-(1/pSum) + 1/pT)*pxSum)/pzSum;//dH/dpx
                matrix(2,5+col*7) =     ((-(1/pSum) + 1/pT)*pySum)/pzSum;//dH/dpy
                matrix(2,6+col*7) =     -(1/pSum) + (pSum - pT)/pow(pzSum,2);//dH/dpz
                //dH/dm=0               

                col++;
        }
        
        return matrix;
}

AlgebraicMatrix MultiTrackPointingKinematicConstraint::positionDerivative(const std::vector<KinematicState> states, const GlobalPoint& point) const{
        int num = states.size();
        if(num<2) throw VertexException("MultiTrackPointingKinematicConstraint::positionDerivative <2 states passed");
        
        //2 equations (for all tracks)
        AlgebraicMatrix  matrix(2,3,0);
        double dx = point.x() - refPoint.x();
        double dy = point.y() - refPoint.y();
        double dz = point.z() - refPoint.z();
        double dT = sqrt(pow(dx,2) + pow(dy,2));
        double ds = sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2));
        
        double pxSum=0, pySum=0, pzSum=0, aSum = 0;
        for(std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++){
                double a = - i->particleCharge() * i->magneticField()->inInverseGeV(i->globalPosition()).z();
                aSum += a;
                
                pxSum += i->kinematicParameters()(3) - a*(point.y() - i->kinematicParameters()(1));
                pySum += i->kinematicParameters()(4) + a*(point.x() - i->kinematicParameters()(0));
                pzSum += i->kinematicParameters()(5);
        }
        double pT = sqrt(pow(pxSum,2) + pow(pySum,2));
        double pSum = sqrt(pow(pxSum,2) + pow(pySum,2) + pow(pzSum,2));

        matrix(1,1) = (-1 + dx/dT)/dy - aSum/pT + (aSum*(pT - pxSum))/pow(pySum,2);//dH/dxv
        matrix(1,2) = 1/dT + (-dT + dx)/pow(dy,2) - (aSum*(pT - pxSum))/(pT*pySum);//dH/dyv
        //dH/dzv=0
        matrix(2,1) = ((1/ds - 1/dT)*dx)/dz + (aSum*(pSum - pT)*pySum)/(pSum*pT*pzSum);//dH/dxv
        matrix(2,2) = ((1/ds - 1/dT)*dy)/dz - (aSum*(pSum - pT)*pxSum)/(pSum*pT*pzSum);//dH/dyv
        matrix(2,3) = 1/ds + (-ds + dT)/pow(dz,2);//dH/dzv
        
        return matrix;
}

int MultiTrackPointingKinematicConstraint::numberOfEquations() const{
        return 2;
}