CurvilinearParametersT.icc

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
///////////////////////////////////////////////////////////////////
// CurvilinearParametersT.icc, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
 
// Gaudi
#include "GaudiKernel/MsgStream.h"
// Trk
#include "TrkEventPrimitives/ParamDefs.h"
#include <cmath>
 
namespace Trk {
// Constructor with TP arguments
template<int DIM, class T, class S>
CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
  const AmgVector(DIM + 2) & parameters,
  std::optional<AmgSymMatrix(DIM)> covariance,
  unsigned int cIdentifier)
  : ParametersBase<DIM, T>(std::move(covariance))
 , m_cIdentifier(cIdentifier)
{
  this->m_position = Amg::Vector3D(parameters[x], parameters[y], parameters[z]);
  this->m_momentum = Amg::Vector3D(parameters[3], parameters[4], parameters[5]);
 
  // flip the charge according to qOverP
  if (parameters[6] < 0.) {
    this->m_chargeDef.setCharge(-1);
    // cppcheck-suppress missingReturn; false positive
  }
  // assign the parameters
  this->m_parameters[locX] = 0.;
  this->m_parameters[locY] = 0.;
  // get phi & theta from the momentum vector
  this->m_parameters[phi] = this->momentum().phi();
  this->m_parameters[theta] = this->momentum().theta();
  this->m_parameters[qOverP] = parameters[6] / this->momentum().mag();
 
  /* we need all the above to be there for the surfac*/
  this->m_surface = S(this->m_position, curvilinearFrame());
}
 
// Constructor with TP arguments
template<int DIM, class T, class S>
CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
  const Amg::Vector3D& pos,
  double tphi,
  double ttheta,
  double tqOverP,
  std::optional<AmgSymMatrix(DIM)> cov,
  unsigned int cIdentifier)
  : ParametersBase<DIM, T>(std::move(cov))
  , m_cIdentifier(cIdentifier)
{
  this->m_position = pos;
  // flip the charge according to qOverP
  if (tqOverP < 0.) {
    this->m_chargeDef.setCharge(-1.);
  } else {
    this->m_chargeDef.setCharge(1.);
  }
 
  // assign the parameters
  this->m_parameters[Trk::locX] = 0.;
  this->m_parameters[Trk::locY] = 0.;
  this->m_parameters[Trk::phi] = tphi;
  this->m_parameters[Trk::theta] = ttheta;
  this->m_parameters[Trk::qOverP] = tqOverP;
 
  // make sure that the position & momentum are calculated
  double p = std::abs(1. / tqOverP);
  this->m_momentum = Amg::Vector3D(p * std::cos(tphi) * std::sin(ttheta),
                                   p * std::sin(tphi) * std::sin(ttheta),
                                   p * std::cos(ttheta));
 
  /* we need all the above for the surface*/
  this->m_surface = S(this->m_position, curvilinearFrame());
}
 
// full global constructor
template<int DIM, class T, class S>
CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
  const Amg::Vector3D& pos,
  const Amg::Vector3D& mom,
  double charge,
  std::optional<AmgSymMatrix(DIM)> cov,
  unsigned int cIdentifier)
  : ParametersBase<DIM, T>(std::move(cov))
  , m_cIdentifier(cIdentifier)
{
  this->m_chargeDef.setCharge(charge);
  // assign the parameters
  this->m_parameters[Trk::locX] = 0.;
  this->m_parameters[Trk::locY] = 0.;
  this->m_parameters[Trk::phi] = mom.phi();
  this->m_parameters[Trk::theta] = mom.theta();
 
  if (charge == 0.) {
    charge = 1.; // such that below is 1./mom.mag()
  }
 
  this->m_parameters[Trk::qOverP] = charge / mom.mag();
  this->m_position = pos;
  this->m_momentum = mom;
 
  // we need all the above to create the surface
  m_surface = S(this->m_position, curvilinearFrame());
}
 
/** the curvilinear parameters identifier */
template<int DIM, class T, class S>
inline unsigned int
CurvilinearParametersT<DIM, T, S>::cIdentifier() const
{
  return m_cIdentifier;
}
 
template<int DIM, class T, class S>
inline void
CurvilinearParametersT<DIM, T, S>::setcIdentifier(unsigned int cIdentifier)
{
  m_cIdentifier = cIdentifier;
}
 
/** Test to see if there's a surface there. */
template<int DIM, class T, class S>
inline bool
CurvilinearParametersT<DIM, T, S>::hasSurface() const
{
  return true;
}
 
/** Access to the Surface method */
template<int DIM, class T, class S>
inline const S&
CurvilinearParametersT<DIM, T, S>::associatedSurface() const
{
  return m_surface;
}
 
// equality operator
template<int DIM, class T, class S>
bool
CurvilinearParametersT<DIM, T, S>::operator==(
  const ParametersBase<DIM, T>& rhs) const
{
  // tolerance for comparing matrices
  constexpr double tolerance = 1e-8;
 
  // make sure we compare objects of same type
  if (this->type() != rhs.type()){
    return false;
  }
  const auto pCasted = static_cast<decltype(this)>(&rhs);
  // comparison to myself?
  if (pCasted == this) {
    return true;
  }
 
  // compare identifier
  if (cIdentifier() != pCasted->cIdentifier()) {
    return false;
  }
 
  // compare UVT frame
  CurvilinearUVT local_curvilinearFrame = curvilinearFrame();
  CurvilinearUVT casted_curvilinearFrame = pCasted->curvilinearFrame();
  if (!local_curvilinearFrame.curvU().isApprox(casted_curvilinearFrame.curvU(),
                                               tolerance)) {
    return false;
  }
  if (!local_curvilinearFrame.curvV().isApprox(casted_curvilinearFrame.curvV(),
                                               tolerance)) {
    return false;
  }
  if (!local_curvilinearFrame.curvT().isApprox(casted_curvilinearFrame.curvT(),
                                               tolerance)) {
    return false;
  }
  // compare surfaces
  if (associatedSurface() != pCasted->associatedSurface()) {
    return false;
  }
 
  // return compatibility of base class parts
  return ParametersBase<DIM, T>::operator==(rhs);
}
 
/** clone */
template<int DIM, class T, class S>
CurvilinearParametersT<DIM, T, S>*
CurvilinearParametersT<DIM, T, S>::clone() const
{
  return new CurvilinearParametersT<DIM, T, S>(*this);
}
 
/** Return the ParametersType enum */
template<int DIM, class T, class S>
inline constexpr ParametersType
CurvilinearParametersT<DIM, T, S>::type() const
{
  return Trk::Curvilinear;
}
 
/** Return the Surface Type (check SurfaceType enums)*/
template<int DIM, class T, class S>
inline constexpr SurfaceType
CurvilinearParametersT<DIM, T, S>::surfaceType() const
{
  return S::staticType;
}
 
// Surface return (with on demand construction)
template<int DIM, class T, class S>
Amg::RotationMatrix3D
CurvilinearParametersT<DIM, T, S>::measurementFrame() const
{
  Amg::RotationMatrix3D mFrame;
  // the columnes
  CurvilinearUVT local_curvilinearFrame = curvilinearFrame();
  mFrame.col(0) = local_curvilinearFrame.curvU();
  mFrame.col(1) = local_curvilinearFrame.curvV();
  mFrame.col(2) = local_curvilinearFrame.curvT();
 
  // return the rotation matrix that defines the curvilinear parameters
  return mFrame;
}
 
template<int DIM, class T, class S>
CurvilinearUVT
CurvilinearParametersT<DIM, T, S>::curvilinearFrame() const
{
  CurvilinearUVT curvilinFrame(this->momentum().unit());
  return curvilinFrame;
}
 
// Screen output dump
template<int DIM, class T, class S>
MsgStream&
CurvilinearParametersT<DIM, T, S>::dump(MsgStream& out) const
{
  out << "CurvilinearParametersT parameters:" << std::endl;
  ParametersBase<DIM, T>::dump(out);
 
  return out;
}
 
// Screen output dump
template<int DIM, class T, class S>
std::ostream&
CurvilinearParametersT<DIM, T, S>::dump(std::ostream& out) const
{
  out << "CurvilinearParametersT parameters:" << std::endl;
  ParametersBase<DIM, T>::dump(out);
 
  return out;
}
 
// private updateParametersHelper
template<int DIM, class T, class S>
void
CurvilinearParametersT<DIM, T, S>::updateParametersHelper(
  const AmgVector(DIM) & updatedParameters)
{
  // valid to use != here, because value is either copied or modified,
  bool updateMomentum =
    (updatedParameters[Trk::phi] != this->m_parameters[Trk::phi]) ||
    (updatedParameters[Trk::theta] != this->m_parameters[Trk::theta]) ||
    (updatedParameters[Trk::qOverP] != this->m_parameters[Trk::qOverP]);
 
  // momentum update is needed
  if (updateMomentum) {
    double phi = updatedParameters[Trk::phi];
    double theta = updatedParameters[Trk::theta];
    double p = std::abs(1. / updatedParameters[Trk::qOverP]);
    this->m_chargeDef.setCharge(sgn(updatedParameters[Trk::qOverP]));
    // assign them and update the momentum 3 vector
    this->m_parameters[Trk::phi] = phi;
    this->m_parameters[Trk::theta] = theta;
    this->m_parameters[Trk::qOverP] = updatedParameters[Trk::qOverP];
    this->m_momentum = Amg::Vector3D(p * std::cos(phi) * std::sin(theta),
                                     p * std::sin(phi) * std::sin(theta),
                                     p * std::cos(theta));
  }
 
  // position update if needed -  loc1
  if (updatedParameters[Trk::loc1] != 0.) {
    this->m_position +=
      updatedParameters[Trk::loc1] * curvilinearFrame().curvU();
  }
  // position update if needed -  loc2
  if (updatedParameters[Trk::loc2] != 0.) {
    this->m_position +=
      updatedParameters[Trk::loc2] * curvilinearFrame().curvV();
  }
  // Reset also the surface
  m_surface = S(this->m_position, curvilinearFrame());
}
 
} // end of namespace Trk