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;
}