BasicSingleVertexState.cc

Go to the documentation of this file.

#include "RecoVertex/VertexPrimitives/interface/BasicSingleVertexState.h"
#include "RecoVertex/VertexPrimitives/interface/VertexException.h"
#include <limits>

namespace {
  constexpr double dNaN = std::numeric_limits<double>::quiet_NaN();
}

BasicSingleVertexState::BasicSingleVertexState()
    : thePos(GlobalPoint(0, 0, 0)),
      theTime(dNaN),
      theErr(AlgebraicSymMatrix44()),
      theWeight(AlgebraicSymMatrix44()),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(0.),
      thePosAvailable(false),
      theTimeAvailable(false),
      theErrAvailable(false),
      theWeightAvailable(false),
      theWeightTimesPosAvailable(false),
      valid(false),
      vertexIs4D(false) {}

BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
                                               const GlobalError& posErr,
                                               const double& weightInMix)
    : thePos(pos),
      theTime(dNaN),
      theErr(posErr),
      theWeight(AlgebraicSymMatrix44()),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(weightInMix),
      thePosAvailable(true),
      theTimeAvailable(false),
      theErrAvailable(true),
      theWeightAvailable(false),
      theWeightTimesPosAvailable(false),
      valid(true),
      vertexIs4D(false) {}

BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
                                               const GlobalWeight& posWeight,
                                               const double& weightInMix)
    : thePos(pos),
      theTime(dNaN),
      theErr(AlgebraicSymMatrix44()),
      theWeight(posWeight),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(weightInMix),
      thePosAvailable(true),
      theTimeAvailable(false),
      theErrAvailable(false),
      theWeightAvailable(true),
      theWeightTimesPosAvailable(false),
      valid(true),
      vertexIs4D(false) {}

BasicSingleVertexState::BasicSingleVertexState(const AlgebraicVector3& weightTimesPosition,
                                               const GlobalWeight& posWeight,
                                               const double& weightInMix)
    : thePos(GlobalPoint(0, 0, 0)),
      theTime(dNaN),
      theErr(AlgebraicSymMatrix44()),
      theWeight(posWeight),
      theWeightTimesPos(weightTimesPosition[0], weightTimesPosition[1], weightTimesPosition[2], 0),
      theWeightInMix(weightInMix),
      thePosAvailable(false),
      theTimeAvailable(false),
      theErrAvailable(false),
      theWeightAvailable(true),
      theWeightTimesPosAvailable(true),
      valid(true),
      vertexIs4D(false) {}

// no-offdiags for time
BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
                                               const GlobalError& posErr,
                                               const double time,
                                               const double timeErr,
                                               const double& weightInMix)
    : thePos(pos),
      theTime(time),
      theErr(posErr),
      theWeight(AlgebraicSymMatrix44()),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(weightInMix),
      thePosAvailable(true),
      theTimeAvailable(true),
      theErrAvailable(true),
      theWeightAvailable(false),
      theWeightTimesPosAvailable(false),
      valid(true),
      vertexIs4D(true) {
  // You dumb bastard. It's not a schooner, its a sailboat.
  GlobalError timeErrMat(0., 0., 0., 0., 0., 0., 0., 0., 0., timeErr * timeErr);
  theErr = theErr + timeErrMat;
}

BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
                                               const GlobalWeight& posWeight,
                                               const double time,
                                               const double timeWeight,
                                               const double& weightInMix)
    : thePos(pos),
      theTime(time),
      theErr(AlgebraicSymMatrix44()),
      theWeight(posWeight),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(weightInMix),
      thePosAvailable(true),
      theTimeAvailable(true),
      theErrAvailable(false),
      theWeightAvailable(true),
      theWeightTimesPosAvailable(false),
      valid(true),
      vertexIs4D(true) {
  GlobalWeight timeWeightMat(0., 0., 0., 0., 0., 0., 0., 0., 0., timeWeight);
  theWeight = theWeight + timeWeightMat;
}

BasicSingleVertexState::BasicSingleVertexState(const AlgebraicVector3& weightTimesPosition,
                                               const GlobalWeight& posWeight,
                                               const double weightTimesTime,
                                               const double timeWeight,
                                               const double& weightInMix)
    : thePos(GlobalPoint(0, 0, 0)),
      theTime(dNaN),
      theErr(AlgebraicSymMatrix44()),
      theWeight(posWeight),
      theWeightTimesPos(weightTimesPosition[0], weightTimesPosition[1], weightTimesPosition[2], weightTimesTime),
      theWeightInMix(weightInMix),
      thePosAvailable(false),
      theTimeAvailable(false),
      theErrAvailable(false),
      theWeightAvailable(true),
      theWeightTimesPosAvailable(true),
      valid(true),
      vertexIs4D(true) {
  GlobalWeight timeWeightMat(0., 0., 0., 0., 0., 0., 0., 0., 0., timeWeight);
  theWeight = theWeight + timeWeightMat;
}

// off-diags for time
BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
                                               const double time,
                                               const GlobalError& posTimeErr,  // fully filled 4x4 matrix
                                               const double& weightInMix)
    : thePos(pos),
      theTime(time),
      theErr(posTimeErr),
      theWeight(AlgebraicSymMatrix44()),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(weightInMix),
      thePosAvailable(true),
      theTimeAvailable(true),
      theErrAvailable(true),
      theWeightAvailable(false),
      theWeightTimesPosAvailable(false),
      valid(true),
      vertexIs4D(true) {}

BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
                                               const double time,
                                               const GlobalWeight& posTimeWeight,
                                               const double& weightInMix)
    : thePos(pos),
      theTime(time),
      theErr(AlgebraicSymMatrix44()),
      theWeight(posTimeWeight),
      theWeightTimesPos(AlgebraicVector4()),
      theWeightInMix(weightInMix),
      thePosAvailable(true),
      theTimeAvailable(true),
      theErrAvailable(false),
      theWeightAvailable(true),
      theWeightTimesPosAvailable(false),
      valid(true),
      vertexIs4D(true) {}

BasicSingleVertexState::BasicSingleVertexState(const AlgebraicVector4& weightTimesPosition,
                                               const GlobalWeight& posWeight,
                                               const double& weightInMix)
    : thePos(GlobalPoint(0, 0, 0)),
      theTime(dNaN),
      theErr(AlgebraicSymMatrix44()),
      theWeight(posWeight),
      theWeightTimesPos(weightTimesPosition),
      theWeightInMix(weightInMix),
      thePosAvailable(false),
      theTimeAvailable(false),
      theErrAvailable(false),
      theWeightAvailable(true),
      theWeightTimesPosAvailable(true),
      valid(true),
      vertexIs4D(true) {}

GlobalPoint BasicSingleVertexState::position() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::position::invalid");
  if (!thePosAvailable)
    computePosition();
  return thePos;
}

GlobalError BasicSingleVertexState::error() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::error::invalid");
  if (!theErrAvailable)
    computeError();
  return GlobalError(theErr.matrix());
}

GlobalError BasicSingleVertexState::error4D() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::error4D::invalid");
  if (!theErrAvailable)
    computeError();
  return theErr;
}

double BasicSingleVertexState::time() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::time::invalid");
  if (!theTimeAvailable)
    computePosition();  // time computed with position (4-vector)
  return theTime;
}

double BasicSingleVertexState::timeError() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::timeError::invalid");
  if (!theTimeAvailable)
    computeError();
  return std::sqrt(theErr.matrix4D()(3, 3));
}

GlobalWeight BasicSingleVertexState::weight() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::weight::invalid");
  if (!theWeightAvailable)
    computeWeight();
  return GlobalWeight(theWeight.matrix());
}

GlobalWeight BasicSingleVertexState::weight4D() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::weight4D::invalid");
  if (!theWeightAvailable)
    computeWeight();
  return theWeight;
}

AlgebraicVector3 BasicSingleVertexState::weightTimesPosition() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::weightTimesPosition::invalid");
  if (!theWeightTimesPosAvailable)
    computeWeightTimesPos();
  return AlgebraicVector3(theWeightTimesPos[0], theWeightTimesPos[1], theWeightTimesPos[2]);
}

AlgebraicVector4 BasicSingleVertexState::weightTimesPosition4D() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::weightTimesPosition4D::invalid");
  if (!theWeightTimesPosAvailable)
    computeWeightTimesPos();
  return theWeightTimesPos;
}

double BasicSingleVertexState::weightInMixture() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::weightInMixture::invalid");
  return theWeightInMix;
}
// RefCountedVertexSeed BasicSingleVertexState::seedWithoutTracks() const
// {
//   RefCountedVertexSeed v = new VertexSeed(position(), error());
//   return v;
// }

void BasicSingleVertexState::computePosition() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::computePosition::invalid");
  AlgebraicVector4 pos = error4D().matrix4D() * weightTimesPosition4D();
  thePos = GlobalPoint(pos[0], pos[1], pos[2]);
  theTime = pos[3];
  thePosAvailable = true;
  theTimeAvailable = true;
}

void BasicSingleVertexState::computeError() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::computeError::invalid");
  int ifail;
  if (vertexIs4D) {
    theErr = weight4D().matrix4D().Inverse(ifail);
    if (ifail != 0)
      throw VertexException("BasicSingleVertexState::could not invert weight matrix");
  } else {
    theErr = weight4D().matrix().Inverse(ifail);
    if (ifail != 0)
      throw VertexException("BasicSingleVertexState::could not invert weight matrix");
  }
  theErrAvailable = true;
}

void BasicSingleVertexState::computeWeight() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::computeWeight::invalid");
  int ifail;
  if (vertexIs4D) {
    theWeight = error4D().matrix4D().Inverse(ifail);
    if (ifail != 0)
      throw VertexException("BasicSingleVertexState::could not invert error matrix");
  } else {
    theWeight = error4D().matrix().Inverse(ifail);
    if (ifail != 0)
      throw VertexException("BasicSingleVertexState::could not invert error matrix");
  }
  theWeightAvailable = true;
}

void BasicSingleVertexState::computeWeightTimesPos() const {
  if (!valid)
    throw VertexException("BasicSingleVertexState::computeWeightTimesPos::invalid");
  AlgebraicVector4 pos;
  pos(0) = position().x();
  pos(1) = position().y();
  pos(2) = position().z();
  if (vertexIs4D) {
    pos(3) = theTime;
  } else {
    pos(3) = 0.;
  }
  theWeightTimesPos = weight4D().matrix4D() * pos;
  theWeightTimesPosAvailable = true;
}