#include "TH3.h"
#include "TError.h"

//______________________________________________________________________________
void TH3_Add(TH3 *h, const TH3 *h1, Double_t c1 = 1.0)
{
// Performs the operation: this = this + c1*h1
// if errors are defined (see TH1::Sumw2), errors are also recalculated.
// Note that if h1 has Sumw2 set, Sumw2 is automatically called for this
// if not already set.
//
// IMPORTANT NOTE1: If you intend to use the errors of this histogram later
// you should call Sumw2 before making this operation.
// This is particularly important if you fit the histogram after TH1::Add
//
// IMPORTANT NOTE2: if h1 has a normalisation factor, the normalisation factor
// is used , ie  this = this + c1*factor*h1
// Use the other TH1::Add function if you do not want this feature

   if (!h || !h1) {
      Error("Add","Attempt to add a non-existing histogram");
      return;
   }

   Int_t nbinsx = h->GetNbinsX();
   Int_t nbinsy = h->GetNbinsY();
   Int_t nbinsz = h->GetNbinsZ();
//   - Check histogram compatibility
   if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()) {
      Error("Add","Attempt to add histograms with different number of bins");
      return;
   }
//   - Issue a Warning if histogram limits are different
   if (h->GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
       h->GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
       h->GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
       h->GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
       h->GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
       h->GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
       Warning("Add","Attempt to add histograms with different axis limits");
   }
   if (h->GetDimension() < 2) nbinsy = -1;
   if (h->GetDimension() < 3) nbinsz = -1;

//    Create Sumw2 if h1 has Sumw2 set
   if (h->GetSumw2N() == 0 && h1->GetSumw2N() != 0) h->Sumw2();

//   - Add statistics
   h->SetEntries(c1*h1->GetEntries());

   Stat_t s1[11], s2[11];
   Int_t i;
   for (i=0;i<11;i++) {s1[i] = s2[i] = 0;}
   h->GetStats(s1);
   h1->GetStats(s2);
   for (i=0;i<10;i++) {
      if (i == 1) s1[i] += c1*c1*s2[i];
      else        s1[i] += c1*s2[i];
   }
   h->PutStats(s1);

   h->SetMinimum();
   h->SetMaximum();

//   - Loop on bins (including underflows/overflows)
   Int_t bin, binx, biny, binz;
   Double_t cu;
   Double_t factor =1;
   if (h1->GetNormFactor() != 0) factor = h1->GetNormFactor()/h1->GetSumOfWeights();;
   for (binz=0;binz<=nbinsz+1;binz++) {
      for (biny=0;biny<=nbinsy+1;biny++) {
         for (binx=0;binx<=nbinsx+1;binx++) {
            bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
            cu  = c1*factor*h1->GetBinContent(bin);
            h->AddBinContent(bin,cu);
            if (h->GetSumw2N()) {
               Double_t error1 = factor*h1->GetBinError(bin);
               h->SetBinError(bin,c1*c1*error1*error1);
            }
         }
      }
   }
}