// @(#)root/hist:$Id$
// Author: Axel Naumann (2011-12-13)

/*************************************************************************
 * Copyright (C) 1995-2012, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include "THn.h"

#include "TClass.h"

namespace {
   //______________________________________________________________________________
   //
   // Helper struct to hold one dimension's bin range for THnBinIter.
   //______________________________________________________________________________
   struct CounterRange_t {
      Int_t i;
      Int_t first;
      Int_t last;
      Int_t len;
      Long64_t cellSize;
   };

   //______________________________________________________________________________
   //
   // THnBinIter iterates over all bins of a THn, recursing over all dimensions.
   //______________________________________________________________________________
   class THnBinIter: public ROOT::THnBaseBinIter {
   public:
      THnBinIter(Int_t dim, const TObjArray* axes, const TNDArray* arr,
                 Bool_t respectAxisRange);
      ~THnBinIter() { delete [] fCounter; }

      Long64_t Next(Int_t* coord = 0);
      Int_t GetCoord(Int_t dim) const { return fCounter[dim].i; }
   private:
      THnBinIter(const THnBinIter&); // intentionally unimplemented
      THnBinIter& operator=(const THnBinIter&); // intentionally unimplemented

   public:
      Int_t fNdimensions;
      Long64_t fIndex;
      const TNDArray* fArray;
      CounterRange_t* fCounter;
   };


   //______________________________________________________________________________
   THnBinIter::THnBinIter(Int_t dim, const TObjArray* axes,
                              const TNDArray* arr, Bool_t respectAxisRange):
      ROOT::THnBaseBinIter(respectAxisRange),
      fNdimensions(dim), fIndex(-1), fArray(arr) {
      // Construct a THnBinIter.
      fCounter = new CounterRange_t[dim]();
      for (Int_t i = 0; i < dim; ++i) {
         TAxis *axis = (TAxis*) axes->At(i);
         fCounter[i].len  = axis->GetNbins() + 2;
         fCounter[i].cellSize  = arr->GetCellSize(i);
         if (!respectAxisRange || !axis->TestBit(TAxis::kAxisRange)) {
            fCounter[i].first = 0;
            fCounter[i].last  = fCounter[i].len - 1;
            fCounter[i].i     = 0;
            continue;
         }
         fHaveSkippedBin = kTRUE;
         Int_t min = axis->GetFirst();
         Int_t max = axis->GetLast();
         if (min == 0 && max == 0) {
            // special case where TAxis::SetBit(kAxisRange) and
            // over- and underflow bins are de-selected.
            // first and last are == 0 due to axis12->SetRange(1, axis12->GetNbins());
            min = 1;
            max = axis->GetNbins();
         }
         fCounter[i].first = min;
         fCounter[i].last  = max;
         fCounter[i].i     = min;
         fIndex += fCounter[i].first * fCounter[i].cellSize;
      }
      // First Next() will increment it:
      --fCounter[dim - 1].i;
   }

   //______________________________________________________________________________
   Long64_t THnBinIter::Next(Int_t* coord /*= 0*/) {
      // Return the current linear bin index (in range), then go to the next bin.
      // If all bins have been visited, return -1.
      if (fNdimensions < 0) return -1; // end
      ++fCounter[fNdimensions - 1].i;
      ++fIndex;
      // Wrap around if needed
      for (Int_t d = fNdimensions - 1; d > 0 && fCounter[d].i > fCounter[d].last; --d) {
         // We skip last + 1..size and 0..first - 1, adjust fIndex
         Int_t skippedCells = fCounter[d].len - (fCounter[d].last + 1);
         skippedCells += fCounter[d].first;
         fIndex += skippedCells * fCounter[d].cellSize;
         fCounter[d].i = fCounter[d].first;
         ++fCounter[d - 1].i;
      }
      if (fCounter[0].i > fCounter[0].last) {
         fNdimensions = -1;
         return -1;
      }
      if (coord) {
         for (Int_t d = 0; d < fNdimensions; ++d) {
            coord[d] = fCounter[d].i;
         }
      }
      return fIndex;
   }
} // unnamed namespce



//______________________________________________________________________________
//
//
//    Multidimensional histogram.
//
// Use a THn if you really, really have to store more than three dimensions,
// and if a large fraction of all bins are filled.
// Better alternatives are
//   * THnSparse if a fraction of all bins are filled
//   * TTree
// The major problem of THn is the memory use caused by n-dimensional
// histogramming: a THnD with 8 dimensions and 100 bins per dimension needs
// more than 2.5GB of RAM!
//
// To construct a THn object you must use one of its templated, derived
// classes:
// THnD (typedef for THnT<Double_t>): bin content held by a Double_t,
// THnF (typedef for THnT<Float_t>): bin content held by a Float_t,
// THnL (typedef for THnT<Long_t>): bin content held by a Long_t,
// THnI (typedef for THnT<Int_t>): bin content held by an Int_t,
// THnS (typedef for THnT<Short_t>): bin content held by a Short_t,
// THnC (typedef for THnT<Char_t>): bin content held by a Char_t,
//
// They take name and title, the number of dimensions, and for each dimension
// the number of bins, the minimal, and the maximal value on the dimension's
// axis. A TH2F h("h","h",10, 0., 10., 20, -5., 5.) would correspond to
//   Int_t bins[2] = {10, 20};
//   Double_t xmin[2] = {0., -5.};
//   Double_t xmax[2] = {10., 5.};
//   THnF hn("hn", "hn", 2, bins, min, max);
//
// * Filling
// A THn is filled just like a regular histogram, using
// THn::Fill(x, weight), where x is a n-dimensional Double_t value.
// To take errors into account, Sumw2() must be called before filling the
// histogram.
// Storage is allocated when the first bin content is stored.
//
// * Projections
// The dimensionality of a THn can be reduced by projecting it to
// 1, 2, 3, or n dimensions, which can be represented by a TH1, TH2, TH3, or
// a THn. See the Projection() members. To only project parts of the
// histogram, call
//   hn->GetAxis(12)->SetRange(from_bin, to_bin);
//
// * Conversion from other histogram classes
// The static factory function THn::CreateHn() can be used to create a THn
// from a TH1, TH2, TH3, THnSparse and (for copying) even from a THn. The
// created THn will have compatble storage type, i.e. calling CreateHn() on
// a TH2F will create a THnF.

ClassImp(THn);

//______________________________________________________________________________
THn::THn(const char* name, const char* title,
         Int_t dim, const Int_t* nbins,
         const Double_t* xmin, const Double_t* xmax):
   THnBase(name, title, dim, nbins, xmin, xmax),
   fSumw2(dim, nbins, kTRUE /*overflow*/),
   fCoordBuf() {
   // Construct a THn.
}

//______________________________________________________________________________
THn::~THn()
{
   // Destruct a THn
   delete [] fCoordBuf;
}


//______________________________________________________________________________
ROOT::THnBaseBinIter* THn::CreateIter(Bool_t respectAxisRange) const
{
   // Create an iterator over all bins. Public interface is THnIter.
   return new THnBinIter(GetNdimensions(), GetListOfAxes(), &GetArray(),
                         respectAxisRange);
}

//______________________________________________________________________________
void THn::Sumw2() {
   // Enable calculation of errors
   if (!GetCalculateErrors()) {
      fTsumw2 = 0.;
   }
   // fill sumw2 array with current content
   TNDArray & content = GetArray();
   Long64_t nbins = GetNbins();
   for (Long64_t ibin = 0; ibin < nbins; ++ibin)
      fSumw2.At(ibin) = content.AtAsDouble(ibin);
}


//______________________________________________________________________________
void THn::AllocCoordBuf() const
{
   // Create the coordinate buffer. Outlined to hide allocation
   // from inlined functions.
   fCoordBuf = new Int_t[fNdimensions]();
}

//______________________________________________________________________________
void THn::InitStorage(Int_t* nbins, Int_t /*chunkSize*/)
{
   // Initialize the storage of a histogram created via Init()
   fCoordBuf = new Int_t[fNdimensions]();
   GetArray().Init(fNdimensions, nbins, true /*addOverflow*/);
   fSumw2.Init(fNdimensions, nbins, true /*addOverflow*/);
}

//______________________________________________________________________________
void THn::Reset(Option_t* option /*= ""*/)
{
   // Reset the contents of a THn.
   GetArray().Reset(option);
   fSumw2.Reset(option);
}