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

Go to the documentation of this file.

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


AlgebraicVector  MultiTrackMassKinematicConstraint::value(const std::vector<KinematicState> states,
                        const GlobalPoint& point) const
{
 if(states.size()<nPart) throw VertexException("MultiTrackMassKinematicConstraint::not enough states given");

 double sumEnergy = 0, sumPx=0, sumPy=0., sumPz=0.;

 double a;
 for (unsigned int i=0;i<nPart;++i) {
    a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();

    sumEnergy += states[i].kinematicParameters().energy();
    sumPx += states[i].kinematicParameters()(3) - a*(point.y() - states[i].kinematicParameters()(1));
    sumPy += states[i].kinematicParameters()(4) + a*(point.x() - states[i].kinematicParameters()(0));
    sumPz += states[i].kinematicParameters()(5);
 }

 double j_m = sumPx*sumPx + sumPy*sumPy + sumPz*sumPz;

 AlgebraicVector res(1,0);
 res(1)  = sumEnergy*sumEnergy - j_m - mass*mass;
 return res;
}

AlgebraicMatrix MultiTrackMassKinematicConstraint::parametersDerivative(const std::vector<KinematicState> states,
                                      const GlobalPoint& point) const
{
  if(states.size()<nPart) throw VertexException("MultiTrackMassKinematicConstraint::not enough states given");

  AlgebraicMatrix res(1,states.size()*7,0);

  double sumEnergy = 0, sumPx=0, sumPy=0., sumPz=0.;

 double a;
  for (unsigned int i=0;i<nPart;++i) {
    a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();

    sumEnergy += states[i].kinematicParameters().energy();
    sumPx += states[i].kinematicParameters()(3) - a*(point.y() - states[i].kinematicParameters()(1));
    sumPy += states[i].kinematicParameters()(4) + a*(point.x() - states[i].kinematicParameters()(0));
    sumPz += states[i].kinematicParameters()(5);
  }

  for (unsigned int i=0;i<nPart;++i) {
    a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();

 //x derivatives:
    res(1,1+i*7) = 2*a*sumPy;

 //y derivatives:
    res(1,2+i*7) = -2*a*sumPx;

 //z components:
    res(1,3+i*7)  = 0.;

 //px components:
    res(1,4+i*7)  = 2*sumEnergy/states[i].kinematicParameters().energy()*states[i].kinematicParameters()(3) - 2*sumPx;

 //py components:
    res(1,5+i*7)  = 2*sumEnergy/states[i].kinematicParameters().energy()*states[i].kinematicParameters()(4) - 2*sumPy;

 //pz1 components:
    res(1,6+i*7)  = 2*sumEnergy/states[i].kinematicParameters().energy()*states[i].kinematicParameters()(5) - 2*sumPz;

 //mass components:
    res(1,7+i*7)  = 2*states[i].kinematicParameters().mass()*sumEnergy/states[i].kinematicParameters().energy();
  }
  return res;
}

AlgebraicMatrix MultiTrackMassKinematicConstraint::positionDerivative(const std::vector<KinematicState> states,
                                    const GlobalPoint& point) const
{
  AlgebraicMatrix res(1,3,0);
  if(states.size()<nPart) throw VertexException("MultiTrackMassKinematicConstraint::not enough states given");

  double sumA = 0, sumPx=0, sumPy=0.;

  double a;
  for (unsigned int i=0;i<nPart;++i) {
    a = -states[i].particleCharge() * states[i].magneticField()->inInverseGeV(states[i].globalPosition()).z();
    sumA += a;

    sumPx += states[i].kinematicParameters()(3) - a*(point.y() - states[i].kinematicParameters()(1));
    sumPy += states[i].kinematicParameters()(4) + a*(point.x() - states[i].kinematicParameters()(0));
  }

 //xv component
  res(1,1) = - 2 * sumPy * sumA;

 //yv component
  res(1,2) =   2 * sumPx * sumA;

 //zv component
  res(1,3) = 0.;

  return res;
}