// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file tools.license for terms.

#ifndef tools_histo_c3d
#define tools_histo_c3d

#include "base_cloud"

#include "../mnmx"

#include "h3d"

namespace tools {
namespace histo {

class c3d : public base_cloud {
public:
  static const std::string& s_class() {
    static const std::string s_v("tools::histo::c3d");
    return s_v;
  }
public:
  bool set_title(const std::string&);
  unsigned int dimension() const {return 3;}
  bool reset();
  unsigned int entries() const;
public:
  double sum_of_weights() const;
  bool convert_to_histogram();
  bool is_converted() const;
  bool scale(double);
public:
  bool fill(double,double,double,double = 1);
  double lower_edge_x() const;
  double upper_edge_x() const;
  double lower_edge_y() const;
  double upper_edge_y() const;
  double lower_edge_z() const;
  double upper_edge_z() const;
  double value_x(unsigned int) const;
  double value_y(unsigned int) const;
  double value_z(unsigned int) const;
  double weight(unsigned int) const;
  double mean_x() const;
  double mean_y() const;
  double mean_z() const;
  double rms_x() const;
  double rms_y() const;
  double rms_z() const;
  bool convert(unsigned int,double,double,
               unsigned int,double,double,
               unsigned int,double,double);
  bool convert(const std::vector<double>&,
                       const std::vector<double>&,
                       const std::vector<double>&);
  const histo::h3d& histogram() const;
  bool fill_histogram(histo::h3d& a_histo) const {
    unsigned int number = m_xs.size();
    for(unsigned int index=0;index<number;index++) {
      if(!a_histo.fill(m_xs[index],m_ys[index],m_zs[index],m_ws[index]))
        return false;
    }
    return true;
  }
  bool set_conversion_parameters(unsigned int,double,double,
                                 unsigned int,double,double,
                                 unsigned int,double,double);

  bool set_histogram(h3d* a_histo){ //we take ownership of a_histo.
    reset();
    m_histo = a_histo;
    return true;
  }
public:
  c3d();
  c3d(const std::string&,int=-1);
  virtual ~c3d(){delete m_histo;}
public:
  c3d(const c3d& a_from)
  :base_cloud(a_from)
  ,m_xs(a_from.m_xs)
  ,m_ys(a_from.m_ys)
  ,m_zs(a_from.m_zs)
  ,m_lower_x(a_from.m_lower_x)
  ,m_upper_x(a_from.m_upper_x)
  ,m_lower_y(a_from.m_lower_y)
  ,m_upper_y(a_from.m_upper_y)
  ,m_lower_z(a_from.m_lower_z)
  ,m_upper_z(a_from.m_upper_z)
  ,m_Sxw(a_from.m_Sxw)
  ,m_Sx2w(a_from.m_Sx2w)
  ,m_Syw(a_from.m_Syw)
  ,m_Sy2w(a_from.m_Sy2w)
  ,m_Szw(a_from.m_Szw)
  ,m_Sz2w(a_from.m_Sz2w)
  ,m_cnv_x_num(a_from.m_cnv_x_num)
  ,m_cnv_x_min(a_from.m_cnv_x_min)
  ,m_cnv_x_max(a_from.m_cnv_x_max)
  ,m_cnv_y_num(a_from.m_cnv_y_num)
  ,m_cnv_y_min(a_from.m_cnv_y_min)
  ,m_cnv_y_max(a_from.m_cnv_y_max)
  ,m_cnv_z_num(a_from.m_cnv_z_num)
  ,m_cnv_z_min(a_from.m_cnv_z_min)
  ,m_cnv_z_max(a_from.m_cnv_z_max)
  ,m_histo(0)
  {
    if(a_from.m_histo) {
      m_histo = new histo::h3d(*a_from.m_histo);
    }
  }

  c3d& operator=(const c3d& a_from) {
    base_cloud::operator=(a_from);
    if(&a_from==this) return *this;
    m_xs = a_from.m_xs;
    m_ys = a_from.m_ys;
    m_zs = a_from.m_zs;
    m_lower_x = a_from.m_lower_x;
    m_upper_x = a_from.m_upper_x;
    m_lower_y = a_from.m_lower_y;
    m_upper_y = a_from.m_upper_y;
    m_lower_z = a_from.m_lower_z;
    m_upper_z = a_from.m_upper_z;
    m_Sxw = a_from.m_Sxw;
    m_Sx2w = a_from.m_Sx2w;
    m_Syw = a_from.m_Syw;
    m_Sy2w = a_from.m_Sy2w;
    m_Szw = a_from.m_Szw;
    m_Sz2w = a_from.m_Sz2w;
    m_cnv_x_num = a_from.m_cnv_x_num;
    m_cnv_x_min = a_from.m_cnv_x_min;
    m_cnv_x_max = a_from.m_cnv_x_max;
    m_cnv_y_num = a_from.m_cnv_y_num;
    m_cnv_y_min = a_from.m_cnv_y_min;
    m_cnv_y_max = a_from.m_cnv_y_max;
    m_cnv_z_num = a_from.m_cnv_z_num;
    m_cnv_z_min = a_from.m_cnv_z_min;
    m_cnv_z_max = a_from.m_cnv_z_max;
    delete m_histo;
    m_histo = 0;
    if(a_from.m_histo) {
      m_histo = new histo::h3d(*a_from.m_histo);
    }
    return *this;
  }

protected:
  void clear();
protected:
  std::vector<double> m_xs;
  std::vector<double> m_ys;
  std::vector<double> m_zs;
  double m_lower_x;
  double m_upper_x;
  double m_lower_y;
  double m_upper_y;
  double m_lower_z;
  double m_upper_z;
  double m_Sxw;
  double m_Sx2w;
  double m_Syw;
  double m_Sy2w;
  double m_Szw;
  double m_Sz2w;
  //
  unsigned int m_cnv_x_num;
  double m_cnv_x_min;
  double m_cnv_x_max;
  unsigned int m_cnv_y_num;
  double m_cnv_y_min;
  double m_cnv_y_max;
  unsigned int m_cnv_z_num;
  double m_cnv_z_min;
  double m_cnv_z_max;
  histo::h3d* m_histo;
};

}}

namespace tools {
namespace histo {

inline
c3d::c3d()
:base_cloud(UNLIMITED())
,m_lower_x(0)
,m_upper_x(0)
,m_lower_y(0)
,m_upper_y(0)
,m_lower_z(0)
,m_upper_z(0)
,m_Sxw(0)
,m_Sx2w(0)
,m_Syw(0)
,m_Sy2w(0)
,m_Szw(0)
,m_Sz2w(0)
,m_cnv_x_num(0)
,m_cnv_x_min(0)
,m_cnv_x_max(0)
,m_cnv_y_num(0)
,m_cnv_y_min(0)
,m_cnv_y_max(0)
,m_cnv_z_num(0)
,m_cnv_z_min(0)
,m_cnv_z_max(0)
,m_histo(0)
{}

inline
c3d::c3d(const std::string& a_title,int aLimit)
:base_cloud(aLimit)
,m_lower_x(0)
,m_upper_x(0)
,m_lower_y(0)
,m_upper_y(0)
,m_lower_z(0)
,m_upper_z(0)
,m_Sxw(0)
,m_Sx2w(0)
,m_Syw(0)
,m_Sy2w(0)
,m_Szw(0)
,m_Sz2w(0)
,m_cnv_x_num(0)
,m_cnv_x_min(0)
,m_cnv_x_max(0)
,m_cnv_y_num(0)
,m_cnv_y_min(0)
,m_cnv_y_max(0)
,m_cnv_z_num(0)
,m_cnv_z_min(0)
,m_cnv_z_max(0)
,m_histo(0)
{
  set_title(a_title);
}

inline
bool c3d::is_converted() const {return m_histo ? true : false;}

inline
void c3d::clear(){
  m_lower_x = 0;
  m_upper_x = 0;
  m_lower_y = 0;
  m_upper_y = 0;
  m_lower_z = 0;
  m_upper_z = 0;
  m_Sw = 0;
  m_Sxw = 0;
  m_Sx2w = 0;
  m_Syw = 0;
  m_Sy2w = 0;
  m_Szw = 0;
  m_Sz2w = 0;
  m_xs.clear();
  m_ys.clear();
  m_zs.clear();
  m_ws.clear();
}

inline
bool c3d::convert( 
 unsigned int aBinsX,double aLowerEdgeX,double aUpperEdgeX
,unsigned int aBinsY,double aLowerEdgeY,double aUpperEdgeY 
,unsigned int aBinsZ,double aLowerEdgeZ,double aUpperEdgeZ 
) { 
  if(m_histo) return true; // Done.
  m_histo = new histo::h3d(base_cloud::title(),
                           aBinsX,aLowerEdgeX,aUpperEdgeX,
                           aBinsY,aLowerEdgeY,aUpperEdgeY,
                           aBinsZ,aLowerEdgeZ,aUpperEdgeZ);
  if(!m_histo) return false;
  bool status = fill_histogram(*m_histo);
  clear();
  return status;
}

inline
bool c3d::convert_to_histogram(){
  if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_x_min) ||
      (m_cnv_y_num<=0) || (m_cnv_y_max<=m_cnv_y_min) ||
      (m_cnv_z_num<=0) || (m_cnv_z_max<=m_cnv_z_min) ) {
    double dx = 0.01 * (upper_edge_x() - lower_edge_x())/BINS();
    double dy = 0.01 * (upper_edge_y() - lower_edge_y())/BINS();
    double dz = 0.01 * (upper_edge_z() - lower_edge_z())/BINS();
    return convert(BINS(),lower_edge_x(),upper_edge_x()+dx,
                   BINS(),lower_edge_y(),upper_edge_y()+dy,
                   BINS(),lower_edge_z(),upper_edge_z()+dz);
  } else {
    return convert(m_cnv_x_num,m_cnv_x_min,m_cnv_x_max,
                   m_cnv_y_num,m_cnv_y_min,m_cnv_y_max,
                   m_cnv_z_num,m_cnv_z_min,m_cnv_z_max);
  }
}

inline
bool c3d::set_title(const std::string& a_title){
  m_title = a_title;
  if(m_histo) m_histo->set_title(a_title);
  return true;
}

inline
bool c3d::scale(double a_scale) {
  if(m_histo) {
    return m_histo->scale(a_scale);
  } else {
    unsigned int number = m_ws.size();
    for(unsigned int index=0;index<number;index++) m_ws[index] *= a_scale;
    m_Sw *= a_scale;
    m_Sxw *= a_scale;
    m_Sx2w *= a_scale;
    m_Syw *= a_scale;
    m_Sy2w *= a_scale;
    m_Szw *= a_scale;
    m_Sz2w *= a_scale;
    return true;
  }
}

inline
bool c3d::set_conversion_parameters(
 unsigned int aCnvXnumber,double aCnvXmin,double aCnvXmax
,unsigned int aCnvYnumber,double aCnvYmin,double aCnvYmax
,unsigned int aCnvZnumber,double aCnvZmin,double aCnvZmax
){
  m_cnv_x_num = aCnvXnumber;
  m_cnv_x_min = aCnvXmin;
  m_cnv_x_max = aCnvXmax;
  m_cnv_y_num = aCnvYnumber;
  m_cnv_y_min = aCnvYmin;
  m_cnv_y_max = aCnvYmax;
  m_cnv_z_num = aCnvZnumber;
  m_cnv_z_min = aCnvZmin;
  m_cnv_z_max = aCnvZmax;
  return true;
}


inline
const h3d& c3d::histogram() const { 
  if(!m_histo) const_cast<c3d&>(*this).convert_to_histogram();
  return *m_histo;
}

inline
bool c3d::reset() { 
  clear();
  delete m_histo;
  m_histo = 0;
  return true;
}

inline
bool c3d::fill(double aX,double aY,double aZ,double aW){
  if(!m_histo && (m_limit!=UNLIMITED()) && ((int)m_xs.size()>=m_limit)){
    convert_to_histogram();
  }

  if(m_histo) {
    return m_histo->fill(aX,aY,aZ,aW);
  } else {
    if(m_xs.size()) {
      m_lower_x = mn<double>(aX,m_lower_x);
      m_upper_x = mx<double>(aX,m_upper_x);
    } else {
      m_lower_x = aX;
      m_upper_x = aX;
    }
    if(m_ys.size()) { 
      m_lower_y = mn<double>(aY,m_lower_y);
      m_upper_y = mx<double>(aY,m_upper_y);
    } else {
      m_lower_y = aY;
      m_upper_y = aY;
    }
    if(m_zs.size()) { 
      m_lower_z = mn<double>(aZ,m_lower_z);
      m_upper_z = mx<double>(aZ,m_upper_z);
    } else {
      m_lower_z = aZ;
      m_upper_z = aZ;
    }
    m_xs.push_back(aX);
    m_ys.push_back(aY);
    m_zs.push_back(aZ);
    m_ws.push_back(aW);
    m_Sw += aW;
    double xw = aX * aW;
    m_Sxw += xw;
    m_Sx2w += aX * xw;
    double yw = aY * aW; 
    m_Syw += yw;
    m_Sy2w += aY * yw;
    double zw = aZ * aW;
    m_Szw += zw;
    m_Sz2w += aZ * zw;
    return true;
  }
}

inline
bool c3d::convert( 
 const std::vector<double>& aEdgesX
,const std::vector<double>& aEdgesY
,const std::vector<double>& aEdgesZ
) { 
  if(m_histo) return true;
  m_histo = new histo::h3d(base_cloud::title(),
                           aEdgesX,aEdgesY,aEdgesZ);
  if(!m_histo) return false;
  bool status = fill_histogram(*m_histo);
  clear();
  return status;
}

inline
double c3d::sum_of_weights() const { 
  return (m_histo ? m_histo->sum_bin_heights() : m_Sw);
}

inline
unsigned int c3d::entries() const { 
  return m_histo ? m_histo->all_entries() : m_ws.size();
}

inline
double c3d::lower_edge_x() const { 
  return m_histo ? m_histo->axis_x().lower_edge() : m_lower_x;
}
inline
double c3d::lower_edge_y() const { 
  return m_histo ? m_histo->axis_y().lower_edge() : m_lower_y;
}
inline
double c3d::lower_edge_z() const { 
  return m_histo ? m_histo->axis_z().lower_edge() : m_lower_z;
}
inline
double c3d::upper_edge_x() const { 
  return m_histo ? m_histo->axis_x().upper_edge() : m_upper_x;
}
inline
double c3d::upper_edge_y() const { 
  return m_histo ? m_histo->axis_y().upper_edge() : m_upper_y;
}
inline
double c3d::upper_edge_z() const { 
  return m_histo ? m_histo->axis_z().upper_edge() : m_upper_z;
}
inline
double c3d::value_x(unsigned int aIndex) const { 
  return m_histo ? 0 : m_xs[aIndex];
}
inline
double c3d::value_y(unsigned int aIndex) const { 
  return m_histo ? 0 : m_ys[aIndex];
}
inline
double c3d::value_z(unsigned int aIndex) const { 
  return m_histo ? 0 : m_zs[aIndex];
}
inline
double c3d::weight(unsigned int aIndex) const { 
  return m_histo ? 0 : m_ws[aIndex];
}
inline
double c3d::mean_x() const {
  return m_histo ? m_histo->mean_x() : (m_Sw?m_Sxw/m_Sw:0);
}
inline
double c3d::mean_y() const {
  return m_histo ? m_histo->mean_y() : (m_Sw?m_Syw/m_Sw:0);
}
inline
double c3d::mean_z() const {
  return m_histo ? m_histo->mean_z() : (m_Sw?m_Szw/m_Sw:0);
}
inline
double c3d::rms_x() const {
  double rms = 0; //FIXME nan.
  if(m_histo) {
    rms = m_histo->rms_x();
  } else {
    if(m_Sw==0) {
    } else {
      double mean = m_Sxw / m_Sw;
      rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - mean * mean));
    }
  }
  return rms;
}
inline
double c3d::rms_y() const {
  double rms = 0; //FIXME nan.
  if(m_histo) {
    rms = m_histo->rms_y();
  } else {
    if(m_Sw==0) {
    } else {
      double mean = m_Syw / m_Sw;
      rms = ::sqrt(::fabs( (m_Sy2w / m_Sw) - mean * mean));
    }
  }
  return rms;
}
inline
double c3d::rms_z() const {
  double rms = 0; //FIXME nan.
  if(m_histo) {
    rms = m_histo->rms_z();
  } else {
    if(m_Sw==0) {
    } else {
      double mean = m_Szw / m_Sw;
      rms = ::sqrt(::fabs( (m_Sz2w / m_Sw) - mean * mean));
    }
  }
  return rms;
}

}}

#endif