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

#ifndef tools_histo_slice
#define tools_histo_slice

#include "h1"
#include "h2"
#include "h3"

namespace tools {
namespace histo {

template <class TC,class TN,class TW,class TH>
inline bool fill_slice_x(const h2<TC,TN,TW,TH>& a_from,
                         int aJbeg,int aJend,
                         h1<TC,TN,TW,TH>& a_to) {

  if(!a_from.dimension()) return false;
  
  typedef typename axis<TC>::bn_t bn_t;

  bn_t jbeg;
  if(!a_from.axis_y().in_range_to_absolute_index(aJbeg,jbeg)) return false;
  bn_t jend;
  if(!a_from.axis_y().in_range_to_absolute_index(aJend,jend)) return false;
  if(jbeg>jend) return false;

  if(a_from.axis_x().bins()!=a_to.axis().bins()) return false;

  histo_data<TC,TN,TW> hdata = a_to.get_histo_data();

  bn_t aoffset,offset,jbin;     
  bn_t yoffset = a_from.axis_y().m_offset;
  
  const std::vector<TN>& af_bin_entries = a_from.bins_entries();
  const std::vector<TW>& af_bin_Sw = a_from.bins_sum_w();
  const std::vector<TW>& af_bin_Sw2 = a_from.bins_sum_w2();
  const std::vector< std::vector<TC> >& af_bin_Sxw = a_from.bins_sum_xw();
  const std::vector< std::vector<TC> >& af_bin_Sx2w = a_from.bins_sum_x2w();

  // Fill also the outflow.
  bn_t abins = hdata.m_axes[0].bins()+2;
  for(bn_t aibin=0;aibin<abins;aibin++) {
    //offset1D = ibin
    aoffset = aibin;
    for(jbin=jbeg;jbin<=jend;jbin++) {
      //offset2D = ibin + jbin * yoffset
      // hdata booked with x then :
      offset = aibin + jbin * yoffset;
      // Bin :
      hdata.m_bin_entries[aoffset] += af_bin_entries[offset];
      hdata.m_bin_Sw[aoffset] += af_bin_Sw[offset];
      hdata.m_bin_Sw2[aoffset] += af_bin_Sw2[offset];
      hdata.m_bin_Sxw[aoffset][0] += af_bin_Sxw[offset][0];
      hdata.m_bin_Sx2w[aoffset][0] += af_bin_Sx2w[offset][0];
    }
  }
  a_to.copy_from_data(hdata);
  a_to.update_fast_getters();
  return true;
}

template <class TC,class TN,class TW,class TH>
inline h1<TC,TN,TW,TH>* slice_x(const h2<TC,TN,TW,TH>& a_from,
                             int aJbeg,int aJend,
                             const std::string& a_title) {
  h1<TC,TN,TW,TH>* slice = new h1<TC,TN,TW,TH>(a_title,
                                   a_from.axis_x().bins(),
                                   a_from.axis_x().lower_edge(),
                                   a_from.axis_x().upper_edge());
  if(!fill_slice_x(a_from,aJbeg,aJend,*slice)) {delete slice;return 0;}
  return slice;
}

template <class TC,class TN,class TW,class TH>
inline h1<TC,TN,TW,TH>* projection_x(const h2<TC,TN,TW,TH>& a_from,const std::string& a_title) {
  return slice_x(a_from,axis<TC>::UNDERFLOW_BIN,axis<TC>::OVERFLOW_BIN,a_title);
}

template <class TC,class TN,class TW,class TH>
inline bool fill_slice_y(const h2<TC,TN,TW,TH>& a_from,
                         int aIbeg,int aIend,
                         h1<TC,TN,TW,TH>& a_to) {

  if(!a_from.dimension()) return false;

  typedef typename axis<TC>::bn_t bn_t;

  bn_t ibeg;
  if(!a_from.axis_x().in_range_to_absolute_index(aIbeg,ibeg)) return false;
  bn_t iend;
  if(!a_from.axis_x().in_range_to_absolute_index(aIend,iend)) return false;
  if(ibeg>iend) return false;

  if(a_from.axis_y().bins()!=a_to.axis().bins()) return false;

  histo_data<TC,TN,TW> hdata = a_to.get_histo_data();

  bn_t aibin,aoffset,offset,ibin;     
  bn_t yoffset = a_from.axis_y().m_offset;

  const std::vector<TN>& af_bin_entries = a_from.bins_entries();
  const std::vector<TW>& af_bin_Sw = a_from.bins_sum_w();
  const std::vector<TW>& af_bin_Sw2 = a_from.bins_sum_w2();
  const std::vector< std::vector<TC> >& af_bin_Sxw = a_from.bins_sum_xw();
  const std::vector< std::vector<TC> >& af_bin_Sx2w = a_from.bins_sum_x2w();

  // Fill also the outflow.
  bn_t abins = hdata.m_axes[0].bins()+2;
  for(aibin=0;aibin<abins;aibin++) {
    //offset1D = ibin
    aoffset = aibin;
    for(ibin=ibeg;ibin<=iend;ibin++) {
      //offset2D = ibin + jbin * yoffset
      // hdata booked with y then :
      offset = ibin + aibin * yoffset;
      // Bin :
      hdata.m_bin_entries[aoffset] += af_bin_entries[offset];
      hdata.m_bin_Sw[aoffset] += af_bin_Sw[offset];
      hdata.m_bin_Sw2[aoffset] += af_bin_Sw2[offset];
      hdata.m_bin_Sxw[aoffset][0] += af_bin_Sxw[offset][1];
      hdata.m_bin_Sx2w[aoffset][0] += af_bin_Sx2w[offset][1];
    }
  }
  a_to.copy_from_data(hdata);
  a_to.update_fast_getters();  
  return true;
}

template <class TC,class TN,class TW,class TH>
inline h1<TC,TN,TW,TH>* slice_y(const h2<TC,TN,TW,TH>& a_from,
                             int aIbeg,int aIend,
                             const std::string& a_title) {
  h1<TC,TN,TW,TH>* slice = new h1<TC,TN,TW,TH>(a_title,
                                   a_from.axis_y().bins(),
                                   a_from.axis_y().lower_edge(),
                                   a_from.axis_y().upper_edge());
  if(!fill_slice_y(a_from,aIbeg,aIend,*slice)) {delete slice;return 0;}
  return slice;
}

template <class TC,class TN,class TW,class TH>
inline h1<TC,TN,TW,TH>* projection_y(const h2<TC,TN,TW,TH>& a_from,const std::string& a_title) {
  return slice_y(a_from,axis<TC>::UNDERFLOW_BIN,axis<TC>::OVERFLOW_BIN,a_title);
}


template <class TC,class TN,class TW,class TH>
inline bool fill_slice_yz(const h3<TC,TN,TW,TH>& a_from,
                          int aIbeg,int aIend,h2<TC,TN,TW,TH>& a_to) {

  if(!a_from.dimension()) return false;

  typedef typename axis<TC>::bn_t bn_t;

  bn_t ibeg;
  if(!a_from.axis_x().in_range_to_absolute_index(aIbeg,ibeg)) return false;
  bn_t iend;
  if(!a_from.axis_x().in_range_to_absolute_index(aIend,iend)) return false;
  if(ibeg>iend) return false;

  if(a_from.axis_y().bins()!=a_to.axis_x().bins()) return false;
  if(a_from.axis_z().bins()!=a_to.axis_y().bins()) return false;

  histo_data<TC,TN,TW> hdata = a_to.get_histo_data();
  
  bn_t aibin,ajbin,aoffset,offset,ibin;     
  
  bn_t ayoffset = hdata.m_axes[1].m_offset;
  bn_t yoffset = a_from.axis_y().m_offset;
  bn_t zoffset = a_from.axis_z().m_offset;
  
  bn_t axbins = hdata.m_axes[0].bins()+2;
  bn_t aybins = hdata.m_axes[1].bins()+2;

  const std::vector<TN>& af_bin_entries = a_from.bins_entries();
  const std::vector<TW>& af_bin_Sw = a_from.bins_sum_w();
  const std::vector<TW>& af_bin_Sw2 = a_from.bins_sum_w2();
  const std::vector< std::vector<TC> >& af_bin_Sxw = a_from.bins_sum_xw();
  const std::vector< std::vector<TC> >& af_bin_Sx2w = a_from.bins_sum_x2w();

  // Fill also the outflow.
  for(aibin=0;aibin<axbins;aibin++) {
    for(ajbin=0;ajbin<aybins;ajbin++) {
      //offset2D = ibin + jbin * m_axes[1].m_offset
      aoffset = aibin + ajbin * ayoffset;
      for(ibin=ibeg;ibin<=iend;ibin++) {
        //offset3D = ibin + jbin * m_axes[1].m_offset + kbin*m_axes[2].m_offset;
        // hdata booked with y-z then :
        offset = ibin + aibin * yoffset + ajbin * zoffset;

        // Bin :
        hdata.m_bin_entries[aoffset] += af_bin_entries[offset];
        hdata.m_bin_Sw[aoffset] += af_bin_Sw[offset];
        hdata.m_bin_Sw2[aoffset] += af_bin_Sw2[offset];
        hdata.m_bin_Sxw[aoffset][0] += af_bin_Sxw[offset][1];
        hdata.m_bin_Sxw[aoffset][1] += af_bin_Sxw[offset][2];
        hdata.m_bin_Sx2w[aoffset][0] += af_bin_Sx2w[offset][1];
        hdata.m_bin_Sx2w[aoffset][1] += af_bin_Sx2w[offset][2];
      }
    }
  }
  a_to.copy_from_data(hdata);
  a_to.update_fast_getters();
  return true;
}

template <class TC,class TN,class TW,class TH>
inline bool fill_slice_xy(const h3<TC,TN,TW,TH>& a_from,
                          int aKbeg,int aKend,h2<TC,TN,TW,TH>& a_to) {
  
  if(!a_from.dimension()) return false;

  typedef typename axis<TC>::bn_t bn_t;

  bn_t kbeg;
  if(!a_from.axis_z().in_range_to_absolute_index(aKbeg,kbeg)) return false;
  bn_t kend;
  if(!a_from.axis_z().in_range_to_absolute_index(aKend,kend)) return false;
  if(kbeg>kend) return false;
  
  if(a_from.axis_x().bins()!=a_to.axis_x().bins()) return false;
  if(a_from.axis_y().bins()!=a_to.axis_y().bins()) return false;

  histo_data<TC,TN,TW> hdata = a_to.get_histo_data();
  
  bn_t kbin;
  bn_t aibin,ajbin,aoffset,offset;     
  
  bn_t ayoffset = hdata.m_axes[1].m_offset;
  bn_t yoffset = a_from.axis_y().m_offset;
  bn_t zoffset = a_from.axis_z().m_offset;
  
  bn_t axbins = hdata.m_axes[0].bins()+2;
  bn_t aybins = hdata.m_axes[1].bins()+2;

  const std::vector<TN>& af_bin_entries = a_from.bins_entries();
  const std::vector<TW>& af_bin_Sw = a_from.bins_sum_w();
  const std::vector<TW>& af_bin_Sw2 = a_from.bins_sum_w2();
  const std::vector< std::vector<TC> >& af_bin_Sxw = a_from.bins_sum_xw();
  const std::vector< std::vector<TC> >& af_bin_Sx2w = a_from.bins_sum_x2w();

  // Fill also the outflow.
  for(aibin=0;aibin<axbins;aibin++) {
    for(ajbin=0;ajbin<aybins;ajbin++) {
      //offset2D = ibin + jbin * m_axes[1].m_offset
      aoffset = aibin + ajbin * ayoffset;
      for(kbin=kbeg;kbin<=kend;kbin++) {
        //offset3D = ibin + jbin * m_axes[1].m_offset + kbin*m_axes[2].m_offset;
        // hdata booked with x-y then :
        offset = aibin + ajbin * yoffset + kbin * zoffset;

        // Bin :
        hdata.m_bin_entries[aoffset] += af_bin_entries[offset];
        hdata.m_bin_Sw[aoffset] += af_bin_Sw[offset];
        hdata.m_bin_Sw2[aoffset] += af_bin_Sw2[offset];
        hdata.m_bin_Sxw[aoffset][0] += af_bin_Sxw[offset][0];
        hdata.m_bin_Sxw[aoffset][1] += af_bin_Sxw[offset][1];
        hdata.m_bin_Sx2w[aoffset][0] += af_bin_Sx2w[offset][0];
        hdata.m_bin_Sx2w[aoffset][1] += af_bin_Sx2w[offset][1];
      }
    }
  }
  a_to.copy_from_data(hdata);
  a_to.update_fast_getters();
  return true;
}

template <class TC,class TN,class TW,class TH>
inline bool fill_slice_xz(const h3<TC,TN,TW,TH>& a_from,
                          int aJbeg,int aJend,h2<TC,TN,TW,TH>& a_to) {
  
  if(!a_from.dimension()) return false;

  typedef typename axis<TC>::bn_t bn_t;

  bn_t jbeg;
  if(!a_from.axis_y().in_range_to_absolute_index(aJbeg,jbeg)) return false;
  bn_t jend;
  if(!a_from.axis_y().in_range_to_absolute_index(aJend,jend)) return false;
  if(jbeg>jend) return false;

  if(a_from.axis_x().bins()!=a_to.axis_x().bins()) return false;
  if(a_from.axis_z().bins()!=a_to.axis_y().bins()) return false;
  
  histo_data<TC,TN,TW> hdata = a_to.get_histo_data();

  bn_t aibin,ajbin,aoffset,offset,jbin;     
  
  bn_t ayoffset = hdata.m_axes[1].m_offset;
  bn_t yoffset = a_from.axis_y().m_offset;
  bn_t zoffset = a_from.axis_z().m_offset;
  
  bn_t axbins = hdata.m_axes[0].bins()+2;
  bn_t aybins = hdata.m_axes[1].bins()+2;

  const std::vector<TN>& af_bin_entries = a_from.bins_entries();
  const std::vector<TW>& af_bin_Sw = a_from.bins_sum_w();
  const std::vector<TW>& af_bin_Sw2 = a_from.bins_sum_w2();
  const std::vector< std::vector<TC> >& af_bin_Sxw = a_from.bins_sum_xw();
  const std::vector< std::vector<TC> >& af_bin_Sx2w = a_from.bins_sum_x2w();

  // Fill also the outflow.
  for(aibin=0;aibin<axbins;aibin++) {
    for(ajbin=0;ajbin<aybins;ajbin++) {
      //offset2D = ibin + jbin * m_axes[1].m_offset
      aoffset = aibin + ajbin * ayoffset;
      for(jbin=jbeg;jbin<=jend;jbin++) {
        //offset3D = ibin + jbin * m_axes[1].m_offset + kbin*m_axes[2].m_offset;
        // hdata booked with x-z then :
        offset = aibin + jbin * yoffset + ajbin * zoffset;
  
        // Bin :
        hdata.m_bin_entries[aoffset] += af_bin_entries[offset];
        hdata.m_bin_Sw[aoffset] += af_bin_Sw[offset];
        hdata.m_bin_Sw2[aoffset] += af_bin_Sw2[offset];
        hdata.m_bin_Sxw[aoffset][0] += af_bin_Sxw[offset][0];
        hdata.m_bin_Sxw[aoffset][1] += af_bin_Sxw[offset][2];
        hdata.m_bin_Sx2w[aoffset][0] += af_bin_Sx2w[offset][0];
        hdata.m_bin_Sx2w[aoffset][1] += af_bin_Sx2w[offset][2];
      }
    }
  }
  
  a_to.copy_from_data(hdata);
  a_to.update_fast_getters();

  return true;
}

template <class TC,class TN,class TW,class TH>
inline h2<TC,TN,TW,TH>* slice_xy(const h3<TC,TN,TW,TH>& a_from,
                             int aKbeg,int aKend,
                             const std::string& a_title) {
  h2<TC,TN,TW,TH>* slice = new h2<TC,TN,TW,TH>(a_title,
                                              a_from.axis_x().bins(),
                                              a_from.axis_x().lower_edge(),
                                              a_from.axis_x().upper_edge(),
                                              a_from.axis_y().bins(),
                                              a_from.axis_y().lower_edge(),
                                              a_from.axis_y().upper_edge());
  if(!fill_slice_xy(a_from,aKbeg,aKend,*slice)) {delete slice;return 0;}
  return slice;
}

template <class TC,class TN,class TW,class TH>
inline h2<TC,TN,TW,TH>* slice_yz(const h3<TC,TN,TW,TH>& a_from,
                             int aIbeg,int aIend,
                             const std::string& a_title) {
  h2<TC,TN,TW,TH>* slice = new h2<TC,TN,TW,TH>(a_title,
                                              a_from.axis_y().bins(),
                                              a_from.axis_y().lower_edge(),
                                              a_from.axis_y().upper_edge(),
                                              a_from.axis_z().bins(),
                                              a_from.axis_z().lower_edge(),
                                              a_from.axis_z().upper_edge());
  if(!fill_slice_yz(a_from,aIbeg,aIend,*slice)) {delete slice;return 0;}
  return slice;
}

template <class TC,class TN,class TW,class TH>
inline h2<TC,TN,TW,TH>* slice_xz(const h3<TC,TN,TW,TH>& a_from,
                             int aJbeg,int aJend,
                             const std::string& a_title) {
  h2<TC,TN,TW,TH>* slice = new h2<TC,TN,TW,TH>(a_title,
                                              a_from.axis_x().bins(),
                                              a_from.axis_x().lower_edge(),
                                              a_from.axis_x().upper_edge(),
                                              a_from.axis_z().bins(),
                                              a_from.axis_z().lower_edge(),
                                              a_from.axis_z().upper_edge());
  if(!fill_slice_xz(a_from,aJbeg,aJend,*slice)) {delete slice;return 0;}
  return slice;
}

}}

#endif