// @(#)root/hist:$Name: v4-04-02f $:$Id: TH4D.cxx,v 1.62.4.1 2005/07/06 11:21:33 rdm Exp $
// Author: Rene Brun   27/10/95

/*************************************************************************
 * Copyright (C) 1995-2000, 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 "TH4D.h"
#include "TRandom.h"
#include <iostream>

using namespace std;

ClassImp(TH4D)

//______________________________________________________________________________
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*
//*-*  The 3-D histogram classes derived from the 1-D histogram classes.
//*-*  all operations are supported (fill, fit).
//*-*  Drawing is currently restricted to one single option.
//*-*  A cloud of points is drawn. The number of points is proportional to
//*-*  cell content.
//*-*
//
//  TH4D a 4-D(3+1) histogram with eight bytes per cell (double)
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

//______________________________________________________________________________
TH4D::TH4D()
    :TH3D(), fExtraDimArray(1), fEXAxisTmp(0)
{
   fSumEntries = 0;
   fIsCopied = kFALSE;
   SetBinsLength(27);
}

//______________________________________________________________________________
TH4D::TH4D(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup
                                     ,Int_t nbinsy,Axis_t ylow,Axis_t yup
                                     ,Int_t nbinsz,Axis_t zlow,Axis_t zup
                                     ,Int_t nbinsu,Axis_t ulow,Axis_t uup)
     :TH3D(name,title,nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup),
      fExtraDimArray(fNcells), fEXAxisTmp(0)
{
//*-*-*-*-*-*-*-*-*Normal constructor for fix bin size 4-D histograms*-*-*-*-*
//*-*              ==================================================

   // start garbage collect...
   TCollection::StartGarbageCollection();

   fSumEntries = 0;
   fNbinsU = nbinsu;
   fMinU   = ulow;
   fMaxU   = uup;
   fIsCopied = kFALSE;
   fEXAxisTmp = new TH1DEX(TString(name) += "_ExtraAxis", 
                           TString(title) += "_ExtraAxis", nbinsu, ulow, uup);


}

//______________________________________________________________________________
TH4D::TH4D(const char *name,const char *title,Int_t nbinsx,const Double_t *xbins
                                           ,Int_t nbinsy,const Double_t *ybins
                                           ,Int_t nbinsz,const Double_t *zbins
                                           ,Int_t nbinsu,const Double_t *ubins)
     :TH3D(name,title,nbinsx,xbins,nbinsy,ybins,nbinsz,zbins), 
      fExtraDimArray(fNcells), fEXAxisTmp(0)
{
//*-*-*-*-*-*-*-*-*Normal constructor for fix bin size 4-D histograms*-*-*-*-*
//*-*              ==================================================

   // start garbage collect...
   TCollection::StartGarbageCollection();

   fSumEntries = 0;
   fNbinsU = nbinsu;
   fMinU   = ubins[0];
   fMaxU   = ubins[nbinsu - 1];
   fIsCopied = kFALSE;
   fEXAxisTmp = new TH1DEX(TString(name) += "_ExtraAxis", 
                           TString(title) += "_ExtraAxis", nbinsu, ubins);
}

//______________________________________________________________________________
TH4D::TH4D(const TH4D &h) : TH3D()
{
   // Copy constructor.
   // The list of functions is not copied. (Use Clone if needed)
   Copy((TObject&)h);
}

//______________________________________________________________________________
TH4D::~TH4D()
{
//
// Clean-up all memory which are allocated inside the class.
// Use garbage collector to delete them in order to avoid double-delete 
// of contents when we use it with TFile.
// TFile::Close() calls TCollection::EmptyGarbageCollection() automatically :(
//

   if (fEXAxisTmp) TCollection::GarbageCollect(fEXAxisTmp);
   if (fExtraDimArray.IsOwner()) {
      for (Int_t i=0; i<fExtraDimArray.GetEntries(); i++) {
         if (fExtraDimArray[i]) TCollection::GarbageCollect(fExtraDimArray[i]);
      }
      // now it must give up ownership...
      fExtraDimArray.SetOwner(kFALSE);
   }
   TCollection::EmptyGarbageCollection();
}

//______________________________________________________________________________
void TH4D::Copy(TObject &obj) const
{
// CAUTION:
// The copy function deletes current extra-dimensions! 
//
   Info("Copy", "Copy is called!");
   if (((TH4D&)obj).fExtraDimArray.IsOwner() == kTRUE) {
      Warning("Copy", "Cannot copy to &obj. Call obj.Reset() first");
      return;
   }
   
   TH3D::Copy(obj);
   ((TH4D&)obj).fExtraDimArray.Delete();
   ((TH4D&)obj).fExtraDimArray.AddAll(&fExtraDimArray);
   ((TH4D&)obj).fExtraDimArray.SetOwner(kFALSE);
   ((TH4D&)obj).fSumEntries = fSumEntries;
   ((TH4D&)obj).fNbinsU = fNbinsU;
   ((TH4D&)obj).fMinU   = fMinU;
   ((TH4D&)obj).fMaxU   = fMaxU;
   ((TH4D&)obj).fIsCopied = kTRUE;
}

//______________________________________________________________________________
TH1DEX* TH4D::GetExtraDimensionPtr(Int_t bin) const
{
//
// bin is calculated by x, y, z, namely, GetBin(x,y,z).
// returns a pointer to TH1DEX assosicated with the bin.
// It returns 0 if extra-dimension for [bin] doesn't exist.
//

   TH1DEX* exdim = 0;
   if (fExtraDimArray.GetLast() >= bin && fExtraDimArray[bin]) {
      exdim = (TH1DEX*)fExtraDimArray[bin];
   }

   return exdim;
}

//______________________________________________________________________________
TH1DEX* TH4D::CreateExtraDimension(Int_t bin)
{
//
// Create new TH1DEX extra-dimension associeted with the bin 
// if no extra-dimension is assigned to the bin.
// Returns a pointer to the created TH1DEX.
//

   if (bin < 0) bin = 0;
   if (bin > fNcells) bin = fNcells;
   TH1DEX* exdim = GetExtraDimensionPtr(bin);
   if (exdim) return exdim;

   if (fIsCopied) {
      Error("CreateExtraDimension", "You are copied. Cannot create extra dimension");
   }

   char buf[1024];
   sprintf(buf, "%s_%d", TNamed::GetName(), bin);

   // create new element of TObjArray
   //Info("CreateExtraDimension", "exdim for bin %d will be created!", bin);

   exdim = new TH1DEX(buf, buf, fNbinsU, fMinU, fMaxU);
   exdim->GetXaxis()->SetTitleOffset(fEXAxisTmp->GetXaxis()->GetTitleOffset());
   exdim->GetXaxis()->SetTitle(fEXAxisTmp->GetXaxis()->GetTitle());
   fExtraDimArray.AddAt(exdim, bin);

   //Info("CreateExtraDimension", "exdim %x is created!", exdim);

   //Should set owner here to avoid memory leak, 
   //This is the temporary solution to avoid double-delete triggerd
   //by TCollection::EmptyGarbageCollection which is called by TFile::close 
   //automatically...
   //
   //fExtraDimArray.SetOwner(kTRUE);

   return exdim;
}

//______________________________________________________________________________
Int_t TH4D::BufferEmpty(Int_t action)
{
// Fill histogram with all entries in the buffer.
// action = -1 histogram is reset and refilled from the buffer (called by THistPainter::Paint)
// action =  0 histogram is filled from the buffer
// action =  1 histogram is filled and buffer is deleted
//             The buffer is automatically deleted when the number of entries
//             in the buffer is greater than the number of entries in the histogram
//
// Modifications from TH3D
// call BufferEmpty for all extra dimensions
//

   Int_t size_1 = TH3D::BufferEmpty(action);

   TH1DEX* exdim;
   for (Int_t bin=0; bin<fNcells; bin++) {
      exdim = GetExtraDimensionPtr(bin); 
      if (exdim) exdim->BufferEmpty(action);
   }

   // it might be bug...
   return size_1;
}

//______________________________________________________________________________
Int_t TH4D::BufferFill(Axis_t x, Axis_t y, Axis_t z, Axis_t u, Stat_t w)
{
// Call BufferFill for 3D + 1D(extra dimemsion)
// It calls TH3D::BufferFill(x,y,z,w), then call TH1DEX::BufferFill(u, w)
// for extra dimemsion which is assigned at fExtraDimArray[FindBin(x,y,z)].
//

   Int_t code = TH3D::BufferFill(x, y, z, w); 
   Int_t bin = TH3D::FindBin(x,y,z); 

   TH1DEX *exdim = CreateExtraDimension(bin);
   exdim->BufferFill(u, w);

   // it might be bug to return fixed value...
   return code;
}

//______________________________________________________________________________
Int_t TH4D::Fill(Axis_t x, Axis_t y, Axis_t z, Axis_t u)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y,z by 1 *-*-*-*-*-*-*-*-*
//*-*                  ====================================
//*-* 
//*-* It calls TH3D::Fill(x,y,z,1), then call TH1DEX::Fill(u, 1)
//*-* for extra dimemsion which is assigned at fExtraDimArray[FindBin(x,y,z)].
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (fBuffer) return BufferFill(x,y,z,u,1);


   // increment number of filled bins (for 4D)
   fSumEntries++;

   // fill TH3D
   Int_t bin = TH3D::Fill(x, y, z, 1);
   //Info("Fill", "bin %d, binx %d, biny %d, binz %d is filled", bin, 
   //      fXaxis.FindBin(x), fYaxis.FindBin(y), fZaxis.FindBin(z));

   TH1DEX *exdim = CreateExtraDimension(bin);
   exdim->Fill(u,1.);

   // it might be bug...
   return bin;
}

//______________________________________________________________________________
Int_t TH4D::Fill(Axis_t x, Axis_t y, Axis_t z, Axis_t u, Stat_t w)
{
//*-*-*-*-*-*-*-*-*-*-*Increment cell defined by x,y,z,u by a weight w*-*-*-*-*
//*-*                  =============================================
//*-* It calls TH3D::Fill(x,y,z,w), then call TH1DEX::Fill(u, w)
//*-* for extra dimemsion which is assigned at fExtraDimArray[GetBin(x,y,z)].
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   if (fBuffer) return BufferFill(x,y,z,u,w);

   // increment number of filled bins (for 4D)
   fSumEntries++;

   // fill TH3D
   Int_t bin = TH3D::Fill(x, y, z, w);
   //Info("Fill", "bin %d, binx %d, biny %d, binz %d is filled", bin, 
   //      fXaxis.FindBin(x), fYaxis.FindBin(y), fZaxis.FindBin(z));

   TH1DEX *exdim = CreateExtraDimension(bin);
   exdim->Fill(u,w);

   // it might be bug...
   return bin;
}

//______________________________________________________________________________
void TH4D::FillRandom(TH1 *h, Int_t ntimes)
{
//*-*-*-*-*-*-*Fill histogram following distribution in histogram h*-*-*-*
//*-*          ====================================================
//*-*
//*-*   The distribution contained in the histogram h (TH4D) is integrated
//*-*   over the channel contents.
//*-*   It is normalized to 1.
//*-*   Getting one random number implies:
//*-*     - Generating a random number between 0 and 1 (say r1)
//*-*     - Look in which bin in the normalized integral r1 corresponds to
//*-*     - Fill histogram channel
//*-*   ntimes random numbers are generated
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-*-*-*-*-*

   if (!h) { Error("FillRandom", "Null histogram"); return; }
   if (fDimension != h->GetDimension()) {
      Error("FillRandom", "Histograms with different dimensions"); return;
   }

   if (h->ComputeIntegral() == 0) return;

   TH4D *h4 = (TH4D*)h;
   Int_t loop;
   Axis_t x,y,z,u;
   for (loop=0;loop<ntimes;loop++) {
      h4->GetRandom4(x,y,z,u);
      Fill(x,y,z,u,1.);
   }
}

//______________________________________________________________________________
Double_t TH4D::ComputeIntegral()
{
//  Compute integral (cumulative sum of bins)
//  The result stored in fIntegral is used by the GetRandom functions.
//  This function is automatically called by GetRandom when the fIntegral
//  array does not exist or when the number of entries in the histogram
//  has changed since the previous call to GetRandom.
//  The resulting integral is normalized to 1
//  
//  Modifications from original function:
//  - calculate 4-dimansional (3+1) cumulative
//  - use Long64_t for bin id
//

   Int_t    bin, binx, biny, binz, binu;
   Long64_t ibin;


   // delete previously computed integral (if any)
   if (fIntegral) delete [] fIntegral;

   //   - Allocate space to store the integral and compute integral
   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   Int_t nbinsu = fNbinsU;
   Long64_t nbins  = nbinsx*nbinsy*nbinsz*nbinsu;

   TH1DEX*    exdim = 0;
   Double_t content = 0;

   fIntegral = new Double_t[nbins+2];
   if (nbins+2 > 113404926) {
      Warning("ComputeIntegral", "array size over size(double) * 113404926 would be too large..");
   }
   ibin = 0;
   fIntegral[ibin] = 0;
   for (binz=1;binz<=nbinsz;binz++) {
      for (biny=1;biny<=nbinsy;biny++) {
         for (binx=1;binx<=nbinsx;binx++) {
            bin  = GetBin(binx, biny, binz);
            exdim = GetExtraDimensionPtr(bin);
            if (exdim) {
               for (binu=1; binu<=nbinsu; binu++) {
                  ibin++;  // start from 1
                  content = exdim->GetBinContent(binu);
                  fIntegral[ibin] = fIntegral[ibin-1] + content;
                  //Info("ComputeIntegral",
                  //     "bin %d, binx %d, biny %d, binz %d, binu %d, ibin %lld integral %f", 
                  //     bin, binx, biny, binz, binu, ibin, fIntegral[ibin]);
               }
            } else {
               for (binu=1; binu<=nbinsu; binu++) {
                  ibin++;
                  fIntegral[ibin] = fIntegral[ibin-1]; // no gain
               }
            }
         }
      }
   }

//   - Normalize integral to 1
   if (fIntegral[nbins] == 0 ) {
      Error("ComputeIntegral", "Integral = zero"); return 0;
   }
   for (ibin=1;ibin<=nbins;ibin++)  fIntegral[ibin] /= fIntegral[nbins];
   fIntegral[nbins+1] = fSumEntries;
   return fIntegral[nbins];
}

//______________________________________________________________________________
 Long64_t TH4D::BinarySearch(Long64_t n, const Double_t *array, Double_t value)
{
// Binary search in an array of n values to locate value.
//
// Array is supposed  to be sorted prior to this call.
// If match is found, function returns position of element.
// If no match found, function gives nearest element smaller than value.
//
// Modified from TMath::BinarySearch:
// Find real edge if several match bins found
//      if (value == array[middle-1]) {
//         if (middle-2 > 0) {
//            while(array[middle-1] == array[middle-2]) middle--;
//         }
//         return middle-1;
//      }
// 

   Long64_t nabove, nbelow, middle;
   nabove = n+1;
   nbelow = 0;
   while(nabove-nbelow > 1) {
      middle = (nabove+nbelow)/2;
      if (value == array[middle-1]) {
         if (middle-2 > 0) {
            while(array[middle-1] == array[middle-2]) middle--;
         }
         return middle-1;
      }
      if (value  < array[middle-1]) nabove = middle;
      else                          nbelow = middle;
   }
   return nbelow-1;
}

//______________________________________________________________________________
void TH4D::DecomposeBin(Long64_t ibin, Int_t &binx, Int_t &biny, 
                                       Int_t &binz, Int_t &binu)
{
   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsu = fNbinsU;
   Long64_t nux    = nbinsu*nbinsx;
   Long64_t nuxy   = nux   *nbinsy;

   
   binz = ibin/nuxy;
   biny = (ibin - binz*nuxy)/nux;
   binx = (ibin - binz*nuxy - biny*nux)/nbinsu;
   binu = ibin - (binz*nuxy + biny*nux + binx*nbinsu);

   // Add 1 to binx, biny, binz. Axes loop in ComputeIntegral start from 1 and end 
   // at nxbins (nybins, nzbins) which requires one increment of above calcuration.
   // CAUTION: do not increment binu, both ibin and binu in ComputeIntegral start from 1!
   binx += 1;
   biny += 1;
   binz += 1;
}

//______________________________________________________________________________
Long64_t TH4D::GetNbins() const
{
   Long64_t nbins  = GetNbinsX() * GetNbinsY() * GetNbinsZ() * fNbinsU;
   return nbins;
}

//______________________________________________________________________________
void TH4D::GetRandom4(Axis_t &x, Axis_t &y, Axis_t &z, Axis_t &u)
{
// return 4 random numbers along axis x , y , z and u distributed according
// the cellcontents of a 3+1-dim histogram

   Long64_t nbins  = GetNbins();

   Double_t integral;
   if (fIntegral) {
      if (fIntegral[nbins+1] != fSumEntries) integral = ComputeIntegral();
      else integral = fIntegral[nbins];
   } else {
      integral = ComputeIntegral();
      if (fIntegral == 0) return;
   }
   if (integral == 0) return;

   Float_t r1 = gRandom->Rndm();
   Long64_t ibin = BinarySearch(nbins,fIntegral,r1);
   Int_t binx, biny, binz, binu;

   DecomposeBin(ibin, binx, biny, binz, binu);  //  decompose ibin into each axes bins

   // choose extra dimension
   Int_t bin = GetBin(binx,biny,binz);
   //Info("GetRandom4", 
   //     "choosed exdim for bin %d, binx %d, biny %d, binz %d, binu %d ibin %lld for random number %f",
   //      bin, binx, biny, binz, binu, ibin, r1);
   TH1DEX *exdim = GetExtraDimensionPtr(bin);

   if (exdim) {
      TAxis *uaxis = exdim->GetXaxis();
      u = uaxis->GetBinLowEdge(binu)
         +uaxis->GetBinWidth(binu)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
      x = fXaxis.GetBinLowEdge(binx) + fXaxis.GetBinWidth(binx)*gRandom->Rndm();
      y = fYaxis.GetBinLowEdge(biny) + fYaxis.GetBinWidth(biny)*gRandom->Rndm();
      z = fZaxis.GetBinLowEdge(binz) + fZaxis.GetBinWidth(binz)*gRandom->Rndm();
   } else {
      Fatal("GetRandom4", "TH4D::GetRandom4: extra dimension doesn't exist!");
   }
}

//______________________________________________________________________________
Stat_t TH4D::Integral(Option_t *option) const
{
//Return integral of bin contents. Only bins in the bins range are considered.
// By default the integral is computed as the sum of bin contents in the range.
// if option "width" is specified, the integral is the sum of
// the bin contents multiplied by the bin width in x, y and in z.

   return Integral(fXaxis.GetFirst(),fXaxis.GetLast(),
                   fYaxis.GetFirst(),fYaxis.GetLast(),
                   fZaxis.GetFirst(),fZaxis.GetLast(),
                   fEXAxisTmp->GetXaxis()->GetFirst(),
                   fEXAxisTmp->GetXaxis()->GetLast(),option);
}

//______________________________________________________________________________
Stat_t TH4D::Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, 
                      Int_t binz1, Int_t binz2, Int_t binu1, Int_t binu2,
                      Option_t *option) const
{
//Return integral of bin contents in range [binx1,binx2],[biny1,biny2],[binz1,binz2]
// [binu1, binu2] for a 4-D histogram
// By default the integral is computed as the sum of bin contents in the range.
// if option "width" is specified, the integral is the sum of
// the bin contents multiplied by the bin width in x, y, z and in u.

   Int_t nbinsx = GetNbinsX();
   Int_t nbinsy = GetNbinsY();
   Int_t nbinsz = GetNbinsZ();
   Int_t nbinsu = fNbinsU;
   if (binx1 < 0) binx1 = 0;
   if (binx2 > nbinsx+1) binx2 = nbinsx+1;
   if (biny1 < 0) biny1 = 0;
   if (biny2 > nbinsy+1) biny2 = nbinsy+1;
   if (binz1 < 0) binz1 = 0;
   if (binz2 > nbinsz+1) binz2 = nbinsz+1;
   if (binu1 < 0) binu1 = 0;
   if (binu2 > nbinsu+1) binu2 = nbinsu+1;
   Stat_t integral = 0;

//*-*- Loop on bins in specified range
   TString opt = option;
   opt.ToLower();
   Bool_t width = kFALSE;
   if (opt.Contains("width")) width = kTRUE;
   Int_t bin, binx, biny, binz, binu;
   TH1DEX *exdim;
   for (binz=binz1;binz<=binz2;binz++) {
      for (biny=biny1;biny<=biny2;biny++) {
         for (binx=binx1;binx<=binx2;binx++) {
            bin  = GetBin(binx, biny, binz);
            exdim = GetExtraDimensionPtr(bin);
            if (exdim) {
               for (binu=1; binu<=nbinsu; binu++) {
                  Stat_t content = exdim->GetBinContent(binu);
                  if (width) {
                     integral += content*fXaxis.GetBinWidth(binx)*fYaxis.GetBinWidth(biny)*fZaxis.GetBinWidth(binz)
                                        *exdim->GetXaxis()->GetBinWidth(binu);
                  } else  {
                     integral += content;
                  }
               }
            }
         }
      }
   }
   return integral;
}

//______________________________________________________________________________
Stat_t TH4D::GetBinContent(Int_t bin, Int_t extrabin) const
{
   if (fBuffer) ((TH4D*)this)->BufferEmpty();
   if (bin < 0) bin = 0;
   if (bin >= fNcells) bin = fNcells-1;
   
   Stat_t content = 0;
   TH1DEX *exdim = GetExtraDimensionPtr(bin);
   if (exdim) content = exdim->GetBinContent(extrabin);

   return content;
}

//______________________________________________________________________________
void TH4D::Reset(Option_t *option)
{
//*-*-*-*-*-*-*-*Reset this histogram: contents, errors, etc*-*-*-*-*-*-*-*
//*-*            ===========================================
//*_* It also orphan the ownership of extra dimensions
//*_*

   TH3D::Reset(option);
   TH1DEX* exdim;
   for (Int_t bin=0; bin<fNcells; bin++) {
      exdim = GetExtraDimensionPtr(bin); 
      if (exdim) exdim->Reset(option);
   }
   fExtraDimArray.SetOwner(kFALSE);
}

//______________________________________________________________________________
void TH4D::SetBinContent(Int_t bin, Int_t extrabin, Stat_t content)
{
// Set bin content
   if (bin < 0) return;
   if (bin >= fNcells) return;

   TH1DEX *exdim = CreateExtraDimension(bin);
   exdim->SetBinContent(extrabin, content);
  
   fSumEntries++;
}

//______________________________________________________________________________
TH4D& TH4D::operator=(const TH4D &h1)
{
   if (this != &h1)  ((TH4D&)h1).Copy(*this);
   return *this;
}