///////////////////////////////////////////////////////////////////////////////////
/// \class RAT::DU::Point3D
///
/// \brief 3D position vector in a specified coordinate system
///
/// \author Jeff Tseng (jeff.tseng@physics.ox.ac.uk)
///
/// REVISION HISTORY:\n
///   9 May 2019: first pull request version
///
///////////////////////////////////////////////////////////////////////////////////
#ifndef __RAT_DU_Point3D_hh__
#define __RAT_DU_Point3D_hh__

#include <TObject.h>
#include <TVector3.h>
#include <TMath.h>
#include <string>
#include <vector>

namespace RAT
{

namespace DU
{

class Point3D : public TObject
{
public:

  //--------------- Defined coordinate systems ----------------------

  /// initialize coordinate system offsets from DB
  static void BeginOfRun();

  /// return number of defined coordinate systems
  ///
  /// @return number of defined coordinate systems
  static size_t GetSystemCount() { return sName.size(); }

  /// query whether a coordinate system has a non-zero z offset
  ///
  /// @return true if specified system has a non-zero z offset
  static Bool_t HasZOffset(const size_t index) { return ((1 << index) & sHasZ) != 0; }

  /// get coordinate system's z offset
  ///
  /// @return z offset relative to the default (PSUP) coordinate system
  static Double_t GetZOffset(const size_t index) { return sValZ.at(index); }

  /// get the coordinate system name
  ///
  /// @return name of the specified coordinate system
  static const std::string GetSystemName(const size_t index) { return sName.at(index); }

  /// get a brief description of the coordinate system
  ///
  /// @return description string for coordinate system
  static const std::string GetSystemDescription(const size_t index) { return sDescription.at(index); }

  /// get the system_id for a given system name
  ///
  /// @return the system id
  static size_t GetSystemId(const std::string name);

  /// ROOT utility to print the defined coordinate systems
  static void DumpSystemNames();

  //--------------- Position vectors ----------------------

  /// Construct a Point3D at the origin with the specific coordinate system
  Point3D(size_t id = 0) : TObject(), fSystem(id), fZOffset() {
    CalculateOffsets(fSystem);
    fPos[0] = fPos[1] = fPos[2] = 0.0;
  }

  /// Copy constructor
  Point3D(const Point3D& p) : TObject(), fSystem(p.fSystem), fZOffset(p.fZOffset) {
    // no need to recalculate offsets, since already done for original Point3D object
    fPos[0] = p.fPos[0];
    fPos[1] = p.fPos[1];
    fPos[2] = p.fPos[2];
  }

  /// Construct a Point3D with the specified coordinate system, but location from old Point3D
  Point3D(size_t id, const Point3D& p) : TObject(), fSystem(id), fZOffset() {
    CalculateOffsets(fSystem);
    SetXYZ(p.fSystem, p.fPos[0], p.fPos[1], p.fPos[2]);
  }

  /// Construct a Point3D with the specified position and coordinate system
  Point3D(size_t id, Double_t x, Double_t y, Double_t z) : TObject(), fSystem(id), fZOffset() {
    CalculateOffsets(fSystem);
    fPos[0] = x;
    fPos[1] = y;
    fPos[2] = z;
  }

  /// Set the x coordinate (in the specified coordinate system), translating if necessary
  /// @return reference to self
  Point3D& SetX(size_t id, Double_t x);

  /// Set the y coordinate (in the specified coordinate system), translating if necessary
  /// @return reference to self
  Point3D& SetY(size_t id, Double_t y);

  /// Set the z coordinate (in the specified coordinate system), translating if necessary
  /// @return reference to self
  Point3D& SetZ(size_t id, Double_t z);

  /// Set the x,y,z coordinates in the specified coordinate system, translating if necessary
  /// @return reference to self
  Point3D& SetXYZ(size_t id, Double_t x, Double_t y, Double_t z);

  /// Explicitly change the coordinate system, translating current position if necessary
  /// @return reference to self
  Point3D& SetCoordinateSystem(size_t id); // only way to change system after ctor!

  /// Set to equivalent point in current coordinate system, translating from argument's system if necessary
  /// @return reference to self
  Point3D& operator=(const Point3D& p);

  /// Return coordinate system identifier
  size_t GetCoordinateSystem() const { return fSystem; }

  /// Return current x coordinate
  Double_t X() const { return fPos[0]; }

  /// Return current y coordinate
  Double_t Y() const { return fPos[1]; }

  /// Return current z coordinate
  Double_t Z() const { return fPos[2]; }

  /// Return coordinate (no bounds checking, so should be 0, 1, or 2)
  Double_t operator() (const size_t index) const { return fPos[index]; }

  /// Return coordinate (no bounds checking, so should be 0, 1, or 2)
  Double_t operator[] (const size_t index) const { return fPos[index]; }

  /// Return square of radius of position from origin
  Double_t Mag2() const { return fPos[0]*fPos[0] + fPos[1]*fPos[1] + fPos[2]*fPos[2]; }

  /// Return the radius, relative to origin
  Double_t Mag() const { return TMath::Sqrt(Mag2()); }

  /// Return cosine of the angle relative to the z axis
  Double_t CosTheta() const { return fPos[2] / Mag(); }

  /// Return the angle in radians, relative to the z axis
  Double_t Theta() const { return TMath::ACos(CosTheta()); }

  /// Return the cosine of the azimuthal angle
  Double_t CosPhi() const { return fPos[0] / Mag(); }

  /// Return the sine of the azimuthal angle
  Double_t SinPhi() const { return fPos[1] / Mag(); }

  /// Return the azimuthal angle, in radians
  Double_t Phi() const { return TMath::ATan2(fPos[1], fPos[0]); }

  /// compare whether points are equivalent (translating to same system).
  /// precision default is machine epsilon
  Bool_t IsEquivalentTo(const Point3D& p, Double_t precision = 2.22e-16) const;

  /// Reflect the position through the origin
  /// @return reflected point
  Point3D Reflection() const;

  /// Get a direction vector
  /// @return this - p vector subtraction
  TVector3 GetDirectionFrom(const Point3D& p) const;

  /// Get a new position based on a displacement from current position
  /// @return a new position
  Point3D GetNewPoint(const TVector3& d) const;

  /// Add a displacement to a position
  /// @return reference to self, having added displacement
  Point3D& Move(const TVector3& d);

  /// Add a displacement to a position (synonym for Move())
  /// @return reference to self, having added displacement
  Point3D& operator+=(const TVector3& d) { return Move(d); }

  /// Sum the position vector (for summing to average)
  /// @return reference to self, including the summed Point3D
  Point3D& Accumulate(const Point3D& p);

  /// Sum the position vector (for summing to average), synonym for Accumulate()
  /// @return reference to self, including the summed Point3D
  Point3D& operator+=(const Point3D& p) { return Accumulate(p); }

  /// Scale the position vector down (for averaging)
  /// @return reference to self
  Point3D& Divide(Double_t n);

  /// Scale the position vector down (for averaging), synonym for Divide()
  /// @return reference to self
  Point3D& operator/=(Double_t n) { return Divide(n); }

  ClassDefNV(Point3D, 0);

private:
  // coordinate systems
  static unsigned int sDefined; // bitmap: 1 if system defined
  static unsigned int sHasZ; // bit map: 1 if non-zero z offset
  static std::vector<Double_t> sValZ; // location of origin in default system
  static std::vector<std::string> sName; // name of coordinate system
  static std::vector<std::string> sDescription; // short description

  // mother volume offsets for default coordinate system
  static std::vector<std::string> sDefMothers; // indices of mother volumes
  static std::vector<Double_t> sDefZOffsets; // z offsets of mother volumes

  size_t fSystem;
  std::vector<Double_t> fZOffset;
  Double_t fPos[3];

  void CalculateOffsets(const size_t id) {
    Double_t dz = GetZOffset(id);
    fZOffset.clear();
    for (size_t i = 0; i < GetSystemCount(); ++i) fZOffset.push_back(dz - GetZOffset(i));
  }

}; // class Point3D

} // namespace DU

} // namespace RAT

#endif