bool equals(Double_t n1, Double_t n2, double ERRORLIMIT)
{
   return fabs( n1 - n2 ) > ERRORLIMIT * fabs(n1);
}

bool equals( TH3* h1, TH3* h2,bool debug = false) { 

   double ERRORLIMIT = 1.E-12;
   int differences = 0;
   for ( int i = 0; i <= h1->GetNbinsX() + 1; ++i )
      for ( int j = 0; j <= h1->GetNbinsY() + 1; ++j )
         for ( int h = 0; h <= h1->GetNbinsY() + 1; ++h )
      {
         Double_t x = h1->GetXaxis()->GetBinCenter(i);
         Double_t y = h1->GetYaxis()->GetBinCenter(j);
         Double_t z = h1->GetZaxis()->GetBinCenter(h);

         bool diff = 0; 
         diff +=  equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
         diff +=  equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
         diff += (bool) equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT);
         diff += (bool) equals(z, h2->GetZaxis()->GetBinCenter(h), ERRORLIMIT);
         diff += (bool) equals(h1->GetBinContent(i,j,h), h2->GetBinContent(i,j,h), ERRORLIMIT);
         if (diff>0) differences++;
         
         if (debug)
         {

            
            cout << equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT) << " "
                 << equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT) << " "
                 << equals(z, h2->GetZaxis()->GetBinCenter(h), ERRORLIMIT) << " "
                 << "[" << x << "," << y << "," << z << "]: " 
                 << h1->GetBinContent(i,j,h) << " +/- " << h1->GetBinError(i,j,h) << " | "
                 << h2->GetBinContent(i,j,h) << " +/- " << h2->GetBinError(i,j,h)
                 << " | " << equals(h1->GetBinContent(i,j,h), h2->GetBinContent(i,j,h), ERRORLIMIT)
                 << " "   << equals(h1->GetBinError(i,j,h)  , h2->GetBinError(i,j,h),   ERRORLIMIT);
            if (diff>0) 
               cout << "\t DIFF BIN !!! ";
            cout << std::endl;                  
         }
      }
   
   if (differences>0) {
      cout << "\nFAILED comparison: Histogram have " << differences << " different bins " << std::endl;
      return false; 
   }
   else 
      return true; 
}


void test3DRebin() { 

   // Tests rebin method for 2D Histogram

   Int_t xrebin = 2;
   Int_t yrebin = 3;
   Int_t zrebin = 4;
   double minRange = 1;
   double maxRange = 10;
   TRandom3 r;
   int nEvents = 1000;

   // make the bins of the orginal histo not an exact divider to leave an extra bin
   TH3D* h3d = new TH3D("h3d","Original Histogram", 
                        xrebin * 4, minRange, maxRange, 
                        yrebin * 3, minRange, maxRange,
                        zrebin * 2, minRange, maxRange);

   UInt_t seed = r.GetSeed();
   r.SetSeed(seed);

   for ( Int_t i = 0; i < nEvents; ++i )
      h3d->Fill( r.Uniform( minRange * .9 , maxRange * 1.1 ), 
                 r.Uniform( minRange * .9 , maxRange * 1.1 ),
                 r.Uniform( minRange * .9 , maxRange * 1.1 ));

   TH3D* h3d2 = (TH3D*) h3d->Rebin3D(xrebin,yrebin, zrebin, "h3-rebin");

   // range of rebinned histogram may be different than original one
   TH3D* h3 = new TH3D("test3DRebin", "test3DRebin", 
                       h3d->GetNbinsX() / xrebin, h3d2->GetXaxis()->GetXmin(), h3d2->GetXaxis()->GetXmax(),
                       h3d->GetNbinsY() / yrebin, h3d2->GetYaxis()->GetXmin(), h3d2->GetYaxis()->GetXmax(),
                       h3d->GetNbinsZ() / zrebin, h3d2->GetZaxis()->GetXmin(), h3d2->GetZaxis()->GetXmax() );

   std::cout << r.Rndm() << std::endl;
   for ( Int_t i = 0; i < nEvents; ++i )
      h3->Fill( r.Uniform( minRange * .9 , maxRange * 1.1 ),
                r.Uniform( minRange * .9 , maxRange * 1.1 ), 
                r.Uniform( minRange * .9 , maxRange * 1.1 ));

   //bool ret = equals("TestIntRebin3D", h3d2, h3, cmpOptStats); // | cmpOptDebug);

   bool ret = equals(h3d2, h3, true);
   // delete h3d;
   // delete h3d2;
   // delete h3;
}