BasicVector.h

Go to the documentation of this file.

#ifndef _utl_BasicVector_h_
#define _utl_BasicVector_h_

#include <CLHEP/Geometry/Normal3D.h>
#include <CLHEP/Geometry/Point3D.h>
#include <CLHEP/Geometry/Vector3D.h>
#include <boost/tuple/tuple.hpp>
#include <utl/CoordinateSystem.h>
#include <utl/unconfig.h>
#include <utl/config.h>
#include <utl/GeometryException.h>
#include <utl/ConsecutiveEnumFactory.h>
#include <utl/ThrowPolicy.h>


namespace utl {

  template<class HepVector>
  class BasicVector {

  public:
    typedef HepVector DataType;

    class CoordinateType {
    public:
      virtual ~CoordinateType() = default;
      virtual double X(double p1, double p2, double p3) const = 0;
      virtual double Y(double p1, double p2, double p3) const = 0;
      virtual double Z(double p1, double p2, double p3) const = 0;

      enum Kind {
        eCartesian   = 0,
        eCylindrical = 1,
        eSpherical   = 2
      };
      static const char* const fgKindTags[3];
      typedef
        utl::ConsecutiveEnumFactory<
          Kind,
          eSpherical,
          fgKindTags,
          utl::ThrowPolicy<Kind, CoordinateSystemException>
        > KindCreator;
      static
      const CoordinateType&
      Create(const Kind k)
      {
        switch (k) {
        case eCartesian: return BasicVector<HepVector>::kCartesian;
        case eCylindrical: return BasicVector<HepVector>::kCylindrical;
        case eSpherical: return BasicVector<HepVector>::kSpherical;
        }
        // This shouldn't happen if the switch is kept in sync with the enum;
        // despite that is put here in order to avoid a compiler warning due
        // to the (supposedly) missing return at the end of the function
        // (which is unreacheable under the former condition).
        std::ostringstream e;
        e << "No conversion from enum value " << k << " to utl::BasicVector<HepVector>::CoordinateType";
        throw CoordinateSystemException(e.str());
      }
    };

  protected:

    class CCartesian : public CoordinateType {
    public:
      virtual double X(double x, double /*y*/, double /*z*/) const override { return x; }
      virtual double Y(double /*x*/, double y, double /*z*/) const override { return y; }
      virtual double Z(double /*x*/, double /*y*/, double z) const override { return z; }
    };


    class CCylindrical : public CoordinateType {
    public:
      virtual double X(double rho, double phi, double z) const override;
      virtual double Y(double rho, double phi, double z) const override;
      virtual double Z(double rho, double phi, double z) const override;
    };


    class CSpherical : public CoordinateType {
    public:
      virtual double X(double r, double theta, double phi) const override;
      virtual double Y(double r, double theta, double phi) const override;
      virtual double Z(double r, double theta, double phi) const override;
    };

  public:
    typedef typename boost::tuple<double, double, double> Triple;

  protected:

    BasicVector(const double x, const double y, const double z,
                const CoordinateSystemPtr& coordinateSystem) :
      fCoordinateSystem(coordinateSystem),
      fVector(x, y, z)
    {
      CoordinateSystemValid(coordinateSystem);
    }


    BasicVector(double p1, double p2, double p3,
                const CoordinateSystemPtr& coordinateSystem,
                const CoordinateType& t) :
      fCoordinateSystem(coordinateSystem),
      fVector(t.X(p1, p2, p3), t.Y(p1, p2, p3), t.Z(p1, p2, p3))
    {
      CoordinateSystemValid(coordinateSystem);
    }

  public:
    bool operator==(const BasicVector& v) const
    { v.TransformTo(fCoordinateSystem); return this->fVector == v.fVector; }

    bool operator!=(const BasicVector& v) const
    { return !(*this == v); }

    double GetX(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.x(); }

    double GetY(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.y(); }

    double GetZ(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.z(); }

    Triple
    GetCoordinates(const CoordinateSystemPtr& coordinateSystem)
      const
    {
      TransformTo(coordinateSystem);
      return boost::make_tuple(fVector.x(), fVector.y(), fVector.z());
    }

    CoordinateSystemPtr GetCoordinateSystem() const
    { return fCoordinateSystem; }

    // universally applicable operations

    // Assignment: rely on compiler to construct the correct operator

    BasicVector& operator*=(const double a) { fVector *= a; return *this; }

    BasicVector& operator/=(const double a) { fVector /= a; return *this; }

    double GetTheta(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.theta(); }
    double GetCosTheta(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.cosTheta(); }
    double GetPhi(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.phi(); }
    double GetR(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.mag(); }
    double GetR2(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.mag2(); }
    double GetRho(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.perp(); }
    double GetRho2(const CoordinateSystemPtr& coordinateSystem) const
    { TransformTo(coordinateSystem); return fVector.perp2(); }

    Triple
    GetSphericalCoordinates(const CoordinateSystemPtr& coordinateSystem)
      const
    {
      TransformTo(coordinateSystem);
      return boost::make_tuple(fVector.r(), fVector.theta(), fVector.phi());
    }

    Triple
    GetCylindricalCoordinates(const CoordinateSystemPtr& coordinateSystem)
      const
    {
      TransformTo(coordinateSystem);
      return boost::make_tuple(fVector.perp(), fVector.phi(), fVector.z());
    }


    bool
    IsCloseTo(const BasicVector& v, const double eps = fgEpsilon)
      const
    {
      v.TransformTo(fCoordinateSystem);
      const double d = (fVector - v.fVector).mag2();
      if (d < eps*eps)
        return true;
      else
        return d / (fVector.mag2() + v.fVector.mag2()) < eps*eps;
    }

    void
    TransformTo(const CoordinateSystemPtr& newCoordinateSystem)
      const
    {
      if (!SameCoordinateSystem(newCoordinateSystem)) {
        fVector = TransformedVector(newCoordinateSystem);
        fCoordinateSystem = newCoordinateSystem;
      }
    }

    template<class V>
    bool SameCoordinateSystem(const BasicVector<V>& v) const
    { return fCoordinateSystem == v.fCoordinateSystem; }

    bool SameCoordinateSystem(const CoordinateSystemPtr& coordinateSystem) const
    { return fCoordinateSystem == coordinateSystem; }

    static const CCartesian kCartesian;
    static const CCylindrical kCylindrical;
    static const CSpherical kSpherical;

  protected:
    ~BasicVector() = default;

    BasicVector() :
      fCoordinateSystem(CoordinateSystem::GetRootCoordinateSystem()),
      fVector(0, 0, 0)
    { }

    BasicVector(const DataType& v,
                const CoordinateSystemPtr& coordinateSystem) :
      fCoordinateSystem(coordinateSystem),
      fVector(v)
    { }

    void Normalize() { fVector = fVector.unit(); }

    DataType
    TransformedVector(const CoordinateSystemPtr& newCoordinateSystem) const;

    mutable CoordinateSystemPtr fCoordinateSystem;
    mutable DataType fVector;

    static const double fgEpsilon;

  };


  typedef BasicVector<HepGeom::Normal3D<double>> AxialVectorBase;
  typedef BasicVector<HepGeom::Point3D<double>> PointBase;
  typedef BasicVector<HepGeom::Vector3D<double>> VectorBase;

  extern template class BasicVector<HepGeom::Point3D<double>>;
  extern template class BasicVector<HepGeom::Vector3D<double>>;
  extern template class BasicVector<HepGeom::Normal3D<double>>;

}


#endif