BasicMultiTrajectoryState.cc

Go to the documentation of this file.

#include "TrackingTools/GsfTools/interface/BasicMultiTrajectoryState.h"

#include "DataFormats/GeometrySurface/interface/BoundPlane.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

using namespace SurfaceSideDefinition;

BasicMultiTrajectoryState::BasicMultiTrajectoryState(const std::vector<TSOS>& tsvec)
    : BasicTrajectoryState(tsvec.front().surface()), theStates(tsvec) {
  // only accept planes!!
  const BoundPlane* bp = dynamic_cast<const BoundPlane*>(&tsvec.begin()->surface());
  if UNLIKELY (bp == nullptr)
    throw cms::Exception("LogicError") << "MultiTrajectoryState constructed on cylinder";

  for (auto i = tsvec.begin(); i != tsvec.end(); i++) {
    if UNLIKELY (!i->isValid()) {
      throw cms::Exception("LogicError") << "MultiTrajectoryState constructed with invalid state";
    }
    if UNLIKELY (i->hasError() != tsvec.front().hasError()) {
      throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states with and without errors";
    }
    if UNLIKELY (&i->surface() != &tsvec.front().surface()) {
      throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states with different surfaces";
    }
    if UNLIKELY (i->surfaceSide() != tsvec.front().surfaceSide()) {
      throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states defined before and after material";
    }
    if UNLIKELY (i->localParameters().pzSign() * tsvec.front().localParameters().pzSign() < 0.) {
      throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states with different signs of local p_z";
    }
  }
  //
  combine();
}

void BasicMultiTrajectoryState::rescaleError(double factor) {
  if UNLIKELY (theStates.empty()) {
    edm::LogError("BasicMultiTrajectoryState") << "Trying to rescale errors of empty MultiTrajectoryState!";
    return;
  }

  for (auto& is : theStates)
    is.rescaleError(factor);
  combine();
}

void BasicMultiTrajectoryState::combine() {
  const std::vector<TrajectoryStateOnSurface>& tsos = theStates;

  if UNLIKELY (tsos.empty()) {
    edm::LogError("MultiTrajectoryStateCombiner") << "Trying to collapse empty set of trajectory states!";
    return;
  }

  double pzSign = tsos.front().localParameters().pzSign();
  for (std::vector<TrajectoryStateOnSurface>::const_iterator it = tsos.begin(); it != tsos.end(); it++) {
    if UNLIKELY (it->localParameters().pzSign() != pzSign) {
      edm::LogError("MultiTrajectoryStateCombiner")
          << "Trying to collapse trajectory states with different signs on p_z!";
      return;
    }
  }

  if UNLIKELY (tsos.size() == 1) {
    BasicTrajectoryState::update(tsos.front().weight(),
                                 tsos.front().localParameters(),
                                 tsos.front().localError(),
                                 tsos.front().surface(),
                                 tsos.front().magneticField(),
                                 tsos.front().surfaceSide());
    return;
  }

  double sumw = 0.;
  //int dim = tsos.front().localParameters().vector().num_row();
  AlgebraicVector5 mean;
  AlgebraicSymMatrix55 covarPart1, covarPart2, covtmp;
  for (auto it1 = tsos.begin(); it1 != tsos.end(); it1++) {
    auto weight = it1->weight();
    auto const& param = it1->localParameters().vector();
    sumw += weight;
    mean += weight * param;
    covarPart1 += weight * it1->localError().matrix();
    for (auto it2 = it1 + 1; it2 != tsos.end(); it2++) {
      AlgebraicVector5 diff = param - it2->localParameters().vector();
      ROOT::Math::AssignSym::Evaluate(covtmp, ROOT::Math::TensorProd(diff, diff));
      covarPart2 += (weight * it2->weight()) * covtmp;
    }
  }
  double sumwI = 1.0 / sumw;
  mean *= sumwI;
  covarPart1 *= sumwI;
  covarPart2 *= (sumwI * sumwI);
  AlgebraicSymMatrix55 covar = covarPart1 + covarPart2;

  BasicTrajectoryState::update(sumw,
                               LocalTrajectoryParameters(mean, pzSign),
                               LocalTrajectoryError(covar),
                               tsos.front().surface(),
                               tsos.front().magneticField(),
                               tsos.front().surfaceSide());
}

void BasicMultiTrajectoryState::update(const LocalTrajectoryParameters& p,
                                       const Surface& aSurface,
                                       const MagneticField* field,
                                       const SurfaceSide side) {
  throw cms::Exception("LogicError",
                       "BasicMultiTrajectoryState::update(LocalTrajectoryParameters, Surface, ...) called even if "
                       "canUpdateLocalParameters() is false");
}

void BasicMultiTrajectoryState::update(double weight,
                                       const LocalTrajectoryParameters& p,
                                       const LocalTrajectoryError& err,
                                       const Surface& aSurface,
                                       const MagneticField* field,
                                       const SurfaceSide side) {
  throw cms::Exception("LogicError",
                       "BasicMultiTrajectoryState::update(LocalTrajectoryParameters, LocalTrajectoryError, ...) called "
                       "even if canUpdateLocalParameters() is false");
}