// @(#)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);
}
 THn.cxx:1
 THn.cxx:2
 THn.cxx:3
 THn.cxx:4
 THn.cxx:5
 THn.cxx:6
 THn.cxx:7
 THn.cxx:8
 THn.cxx:9
 THn.cxx:10
 THn.cxx:11
 THn.cxx:12
 THn.cxx:13
 THn.cxx:14
 THn.cxx:15
 THn.cxx:16
 THn.cxx:17
 THn.cxx:18
 THn.cxx:19
 THn.cxx:20
 THn.cxx:21
 THn.cxx:22
 THn.cxx:23
 THn.cxx:24
 THn.cxx:25
 THn.cxx:26
 THn.cxx:27
 THn.cxx:28
 THn.cxx:29
 THn.cxx:30
 THn.cxx:31
 THn.cxx:32
 THn.cxx:33
 THn.cxx:34
 THn.cxx:35
 THn.cxx:36
 THn.cxx:37
 THn.cxx:38
 THn.cxx:39
 THn.cxx:40
 THn.cxx:41
 THn.cxx:42
 THn.cxx:43
 THn.cxx:44
 THn.cxx:45
 THn.cxx:46
 THn.cxx:47
 THn.cxx:48
 THn.cxx:49
 THn.cxx:50
 THn.cxx:51
 THn.cxx:52
 THn.cxx:53
 THn.cxx:54
 THn.cxx:55
 THn.cxx:56
 THn.cxx:57
 THn.cxx:58
 THn.cxx:59
 THn.cxx:60
 THn.cxx:61
 THn.cxx:62
 THn.cxx:63
 THn.cxx:64
 THn.cxx:65
 THn.cxx:66
 THn.cxx:67
 THn.cxx:68
 THn.cxx:69
 THn.cxx:70
 THn.cxx:71
 THn.cxx:72
 THn.cxx:73
 THn.cxx:74
 THn.cxx:75
 THn.cxx:76
 THn.cxx:77
 THn.cxx:78
 THn.cxx:79
 THn.cxx:80
 THn.cxx:81
 THn.cxx:82
 THn.cxx:83
 THn.cxx:84
 THn.cxx:85
 THn.cxx:86
 THn.cxx:87
 THn.cxx:88
 THn.cxx:89
 THn.cxx:90
 THn.cxx:91
 THn.cxx:92
 THn.cxx:93
 THn.cxx:94
 THn.cxx:95
 THn.cxx:96
 THn.cxx:97
 THn.cxx:98
 THn.cxx:99
 THn.cxx:100
 THn.cxx:101
 THn.cxx:102
 THn.cxx:103
 THn.cxx:104
 THn.cxx:105
 THn.cxx:106
 THn.cxx:107
 THn.cxx:108
 THn.cxx:109
 THn.cxx:110
 THn.cxx:111
 THn.cxx:112
 THn.cxx:113
 THn.cxx:114
 THn.cxx:115
 THn.cxx:116
 THn.cxx:117
 THn.cxx:118
 THn.cxx:119
 THn.cxx:120
 THn.cxx:121
 THn.cxx:122
 THn.cxx:123
 THn.cxx:124
 THn.cxx:125
 THn.cxx:126
 THn.cxx:127
 THn.cxx:128
 THn.cxx:129
 THn.cxx:130
 THn.cxx:131
 THn.cxx:132
 THn.cxx:133
 THn.cxx:134
 THn.cxx:135
 THn.cxx:136
 THn.cxx:137
 THn.cxx:138
 THn.cxx:139
 THn.cxx:140
 THn.cxx:141
 THn.cxx:142
 THn.cxx:143
 THn.cxx:144
 THn.cxx:145
 THn.cxx:146
 THn.cxx:147
 THn.cxx:148
 THn.cxx:149
 THn.cxx:150
 THn.cxx:151
 THn.cxx:152
 THn.cxx:153
 THn.cxx:154
 THn.cxx:155
 THn.cxx:156
 THn.cxx:157
 THn.cxx:158
 THn.cxx:159
 THn.cxx:160
 THn.cxx:161
 THn.cxx:162
 THn.cxx:163
 THn.cxx:164
 THn.cxx:165
 THn.cxx:166
 THn.cxx:167
 THn.cxx:168
 THn.cxx:169
 THn.cxx:170
 THn.cxx:171
 THn.cxx:172
 THn.cxx:173
 THn.cxx:174
 THn.cxx:175
 THn.cxx:176
 THn.cxx:177
 THn.cxx:178
 THn.cxx:179
 THn.cxx:180
 THn.cxx:181
 THn.cxx:182
 THn.cxx:183
 THn.cxx:184
 THn.cxx:185
 THn.cxx:186
 THn.cxx:187
 THn.cxx:188
 THn.cxx:189
 THn.cxx:190
 THn.cxx:191
 THn.cxx:192
 THn.cxx:193
 THn.cxx:194
 THn.cxx:195
 THn.cxx:196
 THn.cxx:197
 THn.cxx:198
 THn.cxx:199
 THn.cxx:200
 THn.cxx:201
 THn.cxx:202
 THn.cxx:203
 THn.cxx:204
 THn.cxx:205
 THn.cxx:206
 THn.cxx:207
 THn.cxx:208
 THn.cxx:209
 THn.cxx:210
 THn.cxx:211
 THn.cxx:212
 THn.cxx:213
 THn.cxx:214
 THn.cxx:215
 THn.cxx:216
 THn.cxx:217
 THn.cxx:218
 THn.cxx:219
 THn.cxx:220
 THn.cxx:221
 THn.cxx:222
 THn.cxx:223
 THn.cxx:224
 THn.cxx:225
 THn.cxx:226
 THn.cxx:227
 THn.cxx:228
 THn.cxx:229
 THn.cxx:230
 THn.cxx:231
 THn.cxx:232
 THn.cxx:233
 THn.cxx:234
 THn.cxx:235
 THn.cxx:236