MultiTrajectoryStateMode.cc

Go to the documentation of this file.

#include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
#include "TrackingTools/GsfTools/interface/GetComponents.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "TrackingTools/GsfTools/interface/MultiGaussianStateTransform.h"
#include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
#include "TrackingTools/GsfTools/interface/GaussianSumUtilities1D.h"

#include <iostream>

namespace multiTrajectoryStateMode {

  bool momentumFromModeCartesian(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
    //
    // clear result vector and check validity of the TSOS
    //
    momentum = GlobalVector(0., 0., 0.);
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // 1D mode computation for px, py and pz
    //
    GetComponents comps(tsos);
    auto const& components = comps();
    auto numb = components.size();
    // vectors of components in x, y and z
    std::vector<SingleGaussianState1D> pxStates;
    pxStates.reserve(numb);
    std::vector<SingleGaussianState1D> pyStates;
    pyStates.reserve(numb);
    std::vector<SingleGaussianState1D> pzStates;
    pzStates.reserve(numb);
    // iteration over components
    for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
      // extraction of parameters and variances
      GlobalVector mom(ic->globalMomentum());
      AlgebraicSymMatrix66 cov(ic->cartesianError().matrix());
      pxStates.push_back(SingleGaussianState1D(mom.x(), cov(3, 3), ic->weight()));
      pyStates.push_back(SingleGaussianState1D(mom.y(), cov(4, 4), ic->weight()));
      pzStates.push_back(SingleGaussianState1D(mom.z(), cov(5, 5), ic->weight()));
    }
    //
    // transformation in 1D multi-states and creation of utility classes
    //
    MultiGaussianState1D pxState(pxStates);
    MultiGaussianState1D pyState(pyStates);
    MultiGaussianState1D pzState(pzStates);
    GaussianSumUtilities1D pxUtils(pxState);
    GaussianSumUtilities1D pyUtils(pyState);
    GaussianSumUtilities1D pzUtils(pzState);
    //
    // cartesian momentum vector from modes
    //
    momentum = GlobalVector(pxUtils.mode().mean(), pyUtils.mode().mean(), pzUtils.mode().mean());
    return true;
  }

  bool positionFromModeCartesian(TrajectoryStateOnSurface const& tsos, GlobalPoint& position) {
    //
    // clear result vector and check validity of the TSOS
    //
    position = GlobalPoint(0., 0., 0.);
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // 1D mode computation for x, y and z
    //

    GetComponents comps(tsos);
    auto const& components = comps();
    auto numb = components.size();
    // vectors of components in x, y and z
    std::vector<SingleGaussianState1D> xStates;
    xStates.reserve(numb);
    std::vector<SingleGaussianState1D> yStates;
    yStates.reserve(numb);
    std::vector<SingleGaussianState1D> zStates;
    zStates.reserve(numb);
    // iteration over components
    for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
      // extraction of parameters and variances
      GlobalPoint pos(ic->globalPosition());
      AlgebraicSymMatrix66 cov(ic->cartesianError().matrix());
      xStates.push_back(SingleGaussianState1D(pos.x(), cov(0, 0), ic->weight()));
      yStates.push_back(SingleGaussianState1D(pos.y(), cov(1, 1), ic->weight()));
      zStates.push_back(SingleGaussianState1D(pos.z(), cov(2, 2), ic->weight()));
    }
    //
    // transformation in 1D multi-states and creation of utility classes
    //
    MultiGaussianState1D xState(xStates);
    MultiGaussianState1D yState(yStates);
    MultiGaussianState1D zState(zStates);
    GaussianSumUtilities1D xUtils(xState);
    GaussianSumUtilities1D yUtils(yState);
    GaussianSumUtilities1D zUtils(zState);
    //
    // cartesian position vector from modes
    //
    position = GlobalPoint(xUtils.mode().mean(), yUtils.mode().mean(), zUtils.mode().mean());
    return true;
  }

  bool momentumFromModeLocal(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
    //
    // clear result vector and check validity of the TSOS
    //
    momentum = GlobalVector(0., 0., 0.);
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // mode computation for local co-ordinates q/p, dx/dz, dy/dz
    //
    double qpMode(0);
    double dxdzMode(0);
    double dydzMode(0);
    //
    // first 3 elements of local parameters = q/p, dx/dz, dy/dz
    //
    for (unsigned int iv = 0; iv < 3; ++iv) {
      // extraction of multi-state using helper class
      MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, iv);
      GaussianSumUtilities1D utils(state1D);
      // mode (in case of failure: mean)
      double result = utils.mode().mean();
      if (!utils.modeIsValid())
        result = utils.mean();
      if (iv == 0)
        qpMode = result;
      else if (iv == 1)
        dxdzMode = result;
      else
        dydzMode = result;
    }
    // local momentum vector from dx/dz, dy/dz and q/p + sign of local pz
    LocalVector localP(dxdzMode, dydzMode, 1.);
    localP *= tsos.localParameters().pzSign() / fabs(qpMode) / sqrt(dxdzMode * dxdzMode + dydzMode * dydzMode + 1.);
    // conversion to global coordinates
    momentum = tsos.surface().toGlobal(localP);
    return true;
  }

  bool momentumFromModeQP(TrajectoryStateOnSurface const& tsos, double& momentum) {
    //
    // clear result vector and check validity of the TSOS
    //
    momentum = 0.;
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // mode computation for local co-ordinates q/p, dx/dz, dy/dz
    //
    double qpMode(0);
    //
    // first element of local parameters = q/p
    //
    // extraction of multi-state using helper class
    MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, 0);
    GaussianSumUtilities1D utils(state1D);
    // mode (in case of failure: mean)
    qpMode = utils.mode().mean();
    if (!utils.modeIsValid())
      qpMode = utils.mean();

    momentum = 1. / fabs(qpMode);
    return true;
  }

  bool momentumFromModeP(TrajectoryStateOnSurface const& tsos, double& momentum) {
    //
    // clear result vector and check validity of the TSOS
    //
    momentum = 0.;
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // first element of local parameters = q/p
    //
    // extraction of multi-state using helper class
    MultiGaussianState1D qpMultiState = MultiGaussianStateTransform::multiState1D(tsos, 0);
    std::vector<SingleGaussianState1D> states(qpMultiState.components());
    // transform from q/p to p
    for (unsigned int i = 0; i < states.size(); ++i) {
      SingleGaussianState1D& qpState = states[i];
      double wgt = qpState.weight();
      double qp = qpState.mean();
      double varQp = qpState.variance();
      double p = 1. / fabs(qp);
      double varP = p * p * p * p * varQp;
      states[i] = SingleGaussianState1D(p, varP, wgt);
    }
    MultiGaussianState1D pMultiState(states);
    GaussianSumUtilities1D utils(pMultiState);
    // mode (in case of failure: mean)
    momentum = utils.mode().mean();
    if (!utils.modeIsValid())
      momentum = utils.mean();

    return true;
  }

  bool positionFromModeLocal(TrajectoryStateOnSurface const& tsos, GlobalPoint& position) {
    //
    // clear result vector and check validity of the TSOS
    //
    position = GlobalPoint(0., 0., 0.);
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // mode computation for local co-ordinates x, y
    //
    double xMode(0);
    double yMode(0);
    //
    // last 2 elements of local parameters = x, y
    //
    for (unsigned int iv = 3; iv < 5; ++iv) {
      // extraction of multi-state using helper class
      MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, iv);
      GaussianSumUtilities1D utils(state1D);
      // mode (in case of failure: mean)
      double result = utils.mode().mean();
      if (!utils.modeIsValid())
        result = utils.mean();
      if (iv == 3)
        xMode = result;
      else
        yMode = result;
    }
    // local position vector from x, y
    LocalPoint localP(xMode, yMode, 0.);
    // conversion to global coordinates
    position = tsos.surface().toGlobal(localP);
    return true;
  }

  bool momentumFromModePPhiEta(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
    //
    // clear result vector and check validity of the TSOS
    //
    momentum = GlobalVector(0., 0., 0.);
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return false;
    }
    //
    // 1D mode computation for p, phi, eta
    //
    GetComponents comps(tsos);
    auto const& components = comps();
    auto numb = components.size();
    // vectors of components in p, phi and eta
    std::vector<SingleGaussianState1D> pStates;
    pStates.reserve(numb);
    std::vector<SingleGaussianState1D> phiStates;
    phiStates.reserve(numb);
    std::vector<SingleGaussianState1D> etaStates;
    etaStates.reserve(numb);
    // covariances in cartesian and p-phi-eta and jacobian
    AlgebraicMatrix33 jacobian;
    AlgebraicSymMatrix33 covCart;
    AlgebraicSymMatrix33 covPPhiEta;
    // iteration over components
    for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
      // parameters
      GlobalVector mom(ic->globalMomentum());
      auto px = mom.x();
      auto py = mom.y();
      auto pz = mom.z();
      auto op = 1. / mom.mag();
      auto opt2 = 1. / mom.perp2();
      auto phi = mom.phi();
      auto eta = mom.eta();
      // jacobian
      jacobian(0, 0) = px * op;
      jacobian(0, 1) = py * op;
      jacobian(0, 2) = pz * op;
      jacobian(1, 0) = py * opt2;
      jacobian(1, 1) = -px * opt2;
      jacobian(1, 2) = 0;
      jacobian(2, 0) = px * pz * opt2 * op;
      jacobian(2, 1) = py * pz * opt2 * op;
      jacobian(2, 2) = -op;
      // extraction of the momentum part from the 6x6 cartesian error matrix
      // and conversion to p-phi-eta
      covCart = ic->cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3);
      covPPhiEta = ROOT::Math::Similarity(jacobian, covCart);
      pStates.push_back(SingleGaussianState1D(1 / op, covPPhiEta(0, 0), ic->weight()));
      phiStates.push_back(SingleGaussianState1D(phi, covPPhiEta(1, 1), ic->weight()));
      etaStates.push_back(SingleGaussianState1D(eta, covPPhiEta(2, 2), ic->weight()));
    }
    //
    // transformation in 1D multi-states and creation of utility classes
    //
    MultiGaussianState1D pState(pStates);
    MultiGaussianState1D phiState(phiStates);
    MultiGaussianState1D etaState(etaStates);
    GaussianSumUtilities1D pUtils(pState);
    GaussianSumUtilities1D phiUtils(phiState);
    GaussianSumUtilities1D etaUtils(etaState);
    //
    // parameters from mode (in case of failure: mean)
    //
    auto p = pUtils.modeIsValid() ? pUtils.mode().mean() : pUtils.mean();
    auto phi = phiUtils.modeIsValid() ? phiUtils.mode().mean() : phiUtils.mean();
    auto eta = etaUtils.modeIsValid() ? etaUtils.mode().mean() : etaUtils.mean();
    //   double theta = 2*atan(exp(-eta));
    auto tanth2 = std::exp(-eta);
    auto pt = p * 2 * tanth2 / (1 + tanth2 * tanth2);             // p*sin(theta)
    auto pz = p * (1 - tanth2 * tanth2) / (1 + tanth2 * tanth2);  // p*cos(theta)
    // conversion to a cartesian momentum vector
    momentum = GlobalVector(pt * cos(phi), pt * sin(phi), pz);
    return true;
  }

  int chargeFromMode(TrajectoryStateOnSurface const& tsos) {
    //
    // clear result vector and check validity of the TSOS
    //
    if (!tsos.isValid()) {
      edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
      return 0;
    }
    //
    // mode computation for local co-ordinates q/p
    // extraction of multi-state using helper class
    MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, 0);
    GaussianSumUtilities1D utils(state1D);
    // mode (in case of failure: mean)
    double result = utils.mode().mean();
    if (!utils.modeIsValid())
      result = utils.mean();

    return result > 0. ? 1 : -1;
  }

}  // namespace multiTrajectoryStateMode