From 754dcf93d1bd3c1a9dfc7905a87e3dda5d587144 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Tue, 2 Dec 2008 17:33:33 +0000
Subject: [PATCH] import from branch new version of stresshistogram from David
 - add different bin sizes in X,Y,Z when doing projection - add test for
 THnSparse - add test for Merge and for label axis

git-svn-id: http://root.cern.ch/svn/root/trunk@26597 27541ba8-7e3a-0410-8455-c3a389f83636
---
 test/stressHistogram.cxx | 1330 +++++++++++++++++++++++++++++---------
 1 file changed, 1025 insertions(+), 305 deletions(-)

diff --git a/test/stressHistogram.cxx b/test/stressHistogram.cxx
index 1ec5dc401fc..68f05b0afed 100644
--- a/test/stressHistogram.cxx
+++ b/test/stressHistogram.cxx
@@ -34,22 +34,32 @@ const int nEvents = 1000;
 const int numberOfBins = 10;
 
 enum compareOptions {
-   cmpOptDebug=1,
-   cmpOptNoError=2,
-   cmpOptStats=4
+   cmpOptPrint=1,
+   cmpOptDebug=2,
+   cmpOptNoError=4,
+   cmpOptStats=8
 };
 
+const int defaultEqualOptions = 0; //cmpOptPrint;
+
 TRandom2 r;
 
 typedef bool ( * pointer2Test) ();
 
-// Methods for histogram comparisions (later implemented)
+struct TTestSuite {
+   unsigned int nTests;
+   char suiteName[75];
+   pointer2Test* tests;
+};
 
-int equals(const char* msg, TH1D* h1, TH1D* h2, int options = 0, double ERRORLIMIT = 1E-15);
-int equals(const char* msg, TH2D* h1, TH2D* h2, int options = 0, double ERRORLIMIT = 1E-15);
-int equals(const char* msg, TH3D* h1, TH3D* h2, int options = 0, double ERRORLIMIT = 1E-15);
-int equals(Double_t n1, Double_t n2, double ERRORLIMIT = 1E-15);
-int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT = 1E-15);
+// Methods for histogram comparisions (later implemented)
+void printResult(const char* msg, bool status);
+int equals(const char* msg, TH1D* h1, TH1D* h2, int options = 0, double ERRORLIMIT = 1E-13);
+int equals(const char* msg, TH2D* h1, TH2D* h2, int options = 0, double ERRORLIMIT = 1E-13);
+int equals(const char* msg, TH3D* h1, TH3D* h2, int options = 0, double ERRORLIMIT = 1E-13);
+int equals(const char* msg, THnSparse* h1, THnSparse* h2, int options, double ERRORLIMIT = 1E-13);
+int equals(Double_t n1, Double_t n2, double ERRORLIMIT = 1E-13);
+int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT = 1E-13);
 ostream& operator<<(ostream& out, TH1D* h);
 // old stresHistOpts.cxx file
 
@@ -86,6 +96,39 @@ bool testAdd1()
    return ret;
 }
 
+bool testAddProfile1() 
+{
+   Double_t c1 = r.Rndm();
+   Double_t c2 = r.Rndm();
+
+   TProfile* p1 = new TProfile("t1D1-p1", "p1-Title", numberOfBins, minRange, maxRange);
+   TProfile* p2 = new TProfile("t1D1-p2", "p2-Title", numberOfBins, minRange, maxRange);
+   TProfile* p3 = new TProfile("t1D1-p3", "p3=c1*p1+c2*p2", numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, 1.0);
+      p3->Fill(x, y,  c1);
+   }
+
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, 1.0);
+      p3->Fill(x, y,  c2);
+   }
+
+   TProfile* p4 = new TProfile("t1D1-p4", "p4=c1*p1+p2*c2", numberOfBins, minRange, maxRange);
+   p4->Add(p1, p2, c1, c2);
+
+   bool ret = equals("Add1DProfile1", p3, p4, cmpOptStats, 1E-13);
+   delete p1;
+   delete p2;
+   delete p3;
+   return ret;
+}
+
 bool testAdd2() 
 {
    Double_t c2 = r.Rndm();
@@ -116,6 +159,36 @@ bool testAdd2()
    return ret;
 }
 
+bool testAddProfile2() 
+{
+   Double_t c2 = r.Rndm();
+
+   TProfile* p5 = new TProfile("t1D2-p5", "p5=   p6+c2*p7", numberOfBins, minRange, maxRange);
+   TProfile* p6 = new TProfile("t1D2-p6", "p6-Title", numberOfBins, minRange, maxRange);
+   TProfile* p7 = new TProfile("t1D2-p7", "p7-Title", numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p6->Fill(x, y, 1.0);
+      p5->Fill(x, y, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p7->Fill(x, y, 1.0);
+      p5->Fill(x, y,  c2);
+   }
+
+   p6->Add(p7, c2);
+   
+   bool ret = equals("Add1DProfile2", p5, p6, cmpOptStats, 1E-13);
+   delete p5;
+   delete p7;
+   return ret;
+}
+
 bool testAdd2D1() 
 {
    Double_t c1 = r.Rndm();
@@ -160,6 +233,50 @@ bool testAdd2D1()
    return ret;
 }
 
+bool testAdd2DProfile1() 
+{
+   Double_t c1 = r.Rndm();
+   Double_t c2 = r.Rndm();
+
+   TProfile2D* p1 = new TProfile2D("t2D1-p1", "p1", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   TProfile2D* p2 = new TProfile2D("t2D1-p2", "p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   TProfile2D* p3 = new TProfile2D("t2D1-p3", "p3=c1*p1+c2*p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, z, 1.0);
+      p3->Fill(x, y, z, c1);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, z, 1.0);
+      p3->Fill(x, y, z, c2);
+   }
+
+   TProfile2D* p4 = new TProfile2D("t2D1-p4", "p4=c1*p1+c2*p2", 
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   p4->Add(p1, p2, c1, c2);
+   bool ret = equals("Add2DProfile1", p3, p4, cmpOptStats , 1E-10);
+   delete p1;
+   delete p2;
+   delete p3;
+   return ret;
+}
+
 bool testAdd2D2() 
 {
    Double_t c2 = r.Rndm();
@@ -199,6 +316,45 @@ bool testAdd2D2()
    return ret;
 }
 
+bool testAdd2DProfile2() 
+{
+   Double_t c2 = r.Rndm();
+
+   TProfile2D* p1 = new TProfile2D("t2D2-p1", "p1", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   TProfile2D* p2 = new TProfile2D("t2D2-p2", "p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   TProfile2D* p3 = new TProfile2D("t2D2-p3", "p3=p1+c2*p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, z, 1.0);
+      p3->Fill(x, y, z, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, z, 1.0);
+      p3->Fill(x, y, z,  c2);
+   }
+
+   p1->Add(p2, c2);
+   bool ret = equals("Add2DProfile2", p3, p1, cmpOptStats, 1E-10);
+   delete p2;
+   delete p3;
+   return ret;
+}
+
 bool testAdd3D1() 
 {
    Double_t c1 = r.Rndm();
@@ -249,6 +405,56 @@ bool testAdd3D1()
    return ret;
 }
 
+bool testAdd3DProfile1() 
+{
+   Double_t c1 = r.Rndm();
+   Double_t c2 = r.Rndm();
+
+   TProfile3D* p1 = new TProfile3D("t3D1-p1", "p1", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   TProfile3D* p2 = new TProfile3D("t3D1-p2", "p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   TProfile3D* p3 = new TProfile3D("t3D1-p3", "p3=c1*p1+c2*p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, z, t, 1.0);
+      p3->Fill(x, y, z, t,  c1);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, z, t, 1.0);
+      p3->Fill(x, y, z, t,  c2);
+   }
+
+   TProfile3D* p4 = new TProfile3D("t3D1-p4", "p4=c1*p1+c2*p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   p4->Add(p1, p2, c1, c2);
+   bool ret = equals("Add3DProfile1", p3, p4, cmpOptStats, 1E-10);
+   delete p1;
+   delete p2;
+   delete p3;
+   return ret;
+}
+
 bool testAdd3D2() 
 {
    Double_t c2 = r.Rndm();
@@ -293,6 +499,50 @@ bool testAdd3D2()
    return ret;
 }
 
+bool testAdd3DProfile2() 
+{
+   Double_t c2 = r.Rndm();
+
+   TProfile3D* p1 = new TProfile3D("t3D2-p1", "p1", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   
+   TProfile3D* p2 = new TProfile3D("t3D2-p2", "p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   
+   TProfile3D* p3 = new TProfile3D("t3D2-p3", "p3=p1+c2*p2", 
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, z, t, 1.0);
+      p3->Fill(x, y, z, t, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, z, t, 1.0);
+      p3->Fill(x, y, z, t,  c2);
+   }
+
+   p1->Add(p2, c2);
+   bool ret = equals("Add3DProfile2", p3, p1, cmpOptStats, 1E-10);
+   delete p2;
+   delete p3;
+   return ret;
+}
+
 bool testMul1() 
 {
    Double_t c1 = r.Rndm();
@@ -667,13 +917,10 @@ bool testDivide1()
       h4->SetBinError( bin, sqrt(error) );
    }
 
-//    cout << h2 << endl;
-//    cout << h3 << endl;
-
    return equals("Divide1D1", h1, h4, cmpOptStats/* | cmpOptDebug*/, 1E-13);
 }
 
-bool stressAssign1D()
+bool testAssign1D()
 {
    TH1D* h1 = new TH1D("=1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
 
@@ -692,12 +939,10 @@ bool stressAssign1D()
    return ret;
 }
 
-bool stressAssignProfile1D()
+bool testAssignProfile1D()
 {
    TProfile* p1 = new TProfile("=1D-p1", "p1-Title", numberOfBins, minRange, maxRange);
 
-   //p1->Sumw2();
-
    for ( Int_t e = 0; e < nEvents; ++e ) {
       Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
       Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
@@ -712,7 +957,7 @@ bool stressAssignProfile1D()
    return ret;
 }
 
-bool stressCopyConstructor1D()
+bool testCopyConstructor1D()
 {
    TH1D* h1 = new TH1D("cc1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
 
@@ -730,12 +975,10 @@ bool stressCopyConstructor1D()
    return ret;
 }
 
-bool stressCopyConstructorProfile1D()
+bool testCopyConstructorProfile1D()
 {
    TProfile* p1 = new TProfile("cc1D-p1", "p1-Title", numberOfBins, minRange, maxRange);
 
-   //p1->Sumw2();
-
    for ( Int_t e = 0; e < nEvents; ++e ) {
       Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
       Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
@@ -749,7 +992,7 @@ bool stressCopyConstructorProfile1D()
    return ret;
 }
 
-bool stressClone1D()
+bool testClone1D()
 {
    TH1D* h1 = new TH1D("cl1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
 
@@ -767,12 +1010,10 @@ bool stressClone1D()
    return ret;
 }
 
-bool stressCloneProfile1D()
+bool testCloneProfile1D()
 {
    TProfile* p1 = new TProfile("cl1D-p1", "p1-Title", numberOfBins, minRange, maxRange);
 
-   //p1->Sumw2();
-
    for ( Int_t e = 0; e < nEvents; ++e ) {
       Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
       Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
@@ -786,7 +1027,7 @@ bool stressCloneProfile1D()
    return ret;
 }
 
-bool stressAssign2D()
+bool testAssign2D()
 {
    TH2D* h1 = new TH2D("=2D-h1", "h1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -810,7 +1051,7 @@ bool stressAssign2D()
    return ret;
 }
 
-bool stressAssignProfile2D()
+bool testAssignProfile2D()
 {
    TProfile2D* p1 = new TProfile2D("=2D-p1", "p1-Title", 
                                    numberOfBins, minRange, maxRange,
@@ -834,7 +1075,7 @@ bool stressAssignProfile2D()
 }
 
 
-bool stressCopyConstructor2D()
+bool testCopyConstructor2D()
 {
    TH2D* h1 = new TH2D("cc2D-h1", "h1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -855,7 +1096,7 @@ bool stressCopyConstructor2D()
    return ret;
 }
 
-bool stressCopyConstructorProfile2D()
+bool testCopyConstructorProfile2D()
 {
    TProfile2D* p1 = new TProfile2D("cc2D-p1", "p1-Title", 
                                    numberOfBins, minRange, maxRange,
@@ -875,7 +1116,7 @@ bool stressCopyConstructorProfile2D()
    return ret;
 }
 
-bool stressClone2D()
+bool testClone2D()
 {
    TH2D* h1 = new TH2D("cl2D-h1", "h1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -896,7 +1137,7 @@ bool stressClone2D()
    return ret;
 }
 
-bool stressCloneProfile2D()
+bool testCloneProfile2D()
 {
    TProfile2D* p1 = new TProfile2D("cl2D-p1", "p1-Title", 
                                    numberOfBins, minRange, maxRange,
@@ -916,7 +1157,7 @@ bool stressCloneProfile2D()
    return ret;
 }
 
-bool stressAssign3D()
+bool testAssign3D()
 {
    TH3D* h1 = new TH3D("=3D-h1", "h1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -943,7 +1184,7 @@ bool stressAssign3D()
    return ret;
 }
 
-bool stressAssignProfile3D()
+bool testAssignProfile3D()
 {
    TProfile3D* p1 = new TProfile3D("=3D-p1", "p1-Title", 
                                    numberOfBins, minRange, maxRange,
@@ -969,7 +1210,7 @@ bool stressAssignProfile3D()
    return ret;
 }
 
-bool stressCopyConstructor3D()
+bool testCopyConstructor3D()
 {
    TH3D* h1 = new TH3D("cc3D-h1", "h1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -992,7 +1233,7 @@ bool stressCopyConstructor3D()
    return ret;
 }
 
-bool stressCopyConstructorProfile3D()
+bool testCopyConstructorProfile3D()
 {
    TProfile3D* p1 = new TProfile3D("cc3D-p1", "p1-Title", 
                                    numberOfBins, minRange, maxRange,
@@ -1014,7 +1255,7 @@ bool stressCopyConstructorProfile3D()
    return ret;
 }
 
-bool stressClone3D()
+bool testClone3D()
 {
    TH3D* h1 = new TH3D("cl3D-h1", "h1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -1037,7 +1278,7 @@ bool stressClone3D()
    return ret;
 }
 
-bool stressCloneProfile3D()
+bool testCloneProfile3D()
 {
    TProfile3D* p1 = new TProfile3D("cl3D-p1", "p1-Title", 
                        numberOfBins, minRange, maxRange,
@@ -1086,8 +1327,6 @@ bool testWriteReadProfile1D()
 {
    TProfile* p1 = new TProfile("wr1D-p1", "p1-Title", numberOfBins, minRange, maxRange);
 
-   //p1->Sumw2();
-
    for ( Int_t e = 0; e < nEvents; ++e ) {
       Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
       Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
@@ -1216,51 +1455,332 @@ bool testWriteReadProfile3D()
    return ret;
 }
 
-bool stressHistOpts()
+bool testMerge1D() 
 {
-   r.SetSeed(0);
-   const unsigned int numberOfTests = 36;
-   pointer2Test testPointer[numberOfTests] = {  testAdd1,                testAdd2, 
-                                                testAdd2D1,              testAdd2D2,
-                                                testAdd3D1,              testAdd3D2, 
-                                                testMul1,                testMul2,
-                                                testMul2D1,              testMul2D2,
-                                                testMul3D1,              testMul3D2, 
-                                                //testDivide1,
-                                                stressAssign1D,          stressAssignProfile1D, 
-                                                stressCopyConstructor1D, stressCopyConstructorProfile1D, 
-                                                stressClone1D,           stressCloneProfile1D,
-                                                stressAssign2D,          stressAssignProfile2D,
-                                                stressCopyConstructor2D, stressCopyConstructorProfile2D,
-                                                stressClone2D,           stressCloneProfile2D,
-                                                stressAssign3D,          stressAssignProfile3D,
-                                                stressCopyConstructor3D, stressCopyConstructorProfile3D,
-                                                stressClone3D,           stressCloneProfile3D,
-                                                testWriteRead1D,         testWriteReadProfile1D,
-                                                testWriteRead2D,         testWriteReadProfile2D,
-                                                testWriteRead3D,         testWriteReadProfile3D
-   };
+   TH1D* h1 = new TH1D("merge1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("merge1D-h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("merge1D-h3", "h3-Title", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("merge1D-h4", "h4-Title", numberOfBins, minRange, maxRange);
 
-   // Still to do: testDivide2, testDivide2D1, testDivide2D2 and
-   // testDivide3D1, testDivide3D2. 
+   h1->Sumw2();h2->Sumw2();h3->Sumw2();
 
-   // It depends on whether we can solve the problem with the 1D test
-   // already done. It seems like there is something wrong with the
-   // way the errors are being calculated. We have a formula to
-   // calculate it by hand. Nevertheless the Divide method errors and
-   // the ones calculated by ourselves differ a bit from those (in the
-   // order of 1E-1).
-   
-   bool status = false;
-   for ( unsigned int i = 0; i < numberOfTests; ++i )
-      status |= testPointer[i]();
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h1->Fill(value, 1.0);
+      h4->Fill(value, 1.0);
+   }
 
-   return status;
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h2->Fill(value, 1.0);
+      h4->Fill(value, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h3->Fill(value, 1.0);
+      h4->Fill(value, 1.0);
+   }
+
+   TList *list = new TList;
+   list->Add(h2);
+   list->Add(h3);
+
+   h1->Merge(list);
+
+   bool ret = equals("Merge1D", h1, h4, cmpOptStats, 1E-10);
+   delete h1;
+   delete h2;
+   delete h3;
+   return ret;
+}
+
+bool testMergeProf1D() 
+{
+   TProfile* p1 = new TProfile("merge1D-p1", "p1-Title", numberOfBins, minRange, maxRange);
+   TProfile* p2 = new TProfile("merge1D-p2", "p2-Title", numberOfBins, minRange, maxRange);
+   TProfile* p3 = new TProfile("merge1D-p3", "p3-Title", numberOfBins, minRange, maxRange);
+   TProfile* p4 = new TProfile("merge1D-p4", "p4-Title", numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, 1.0);
+      p4->Fill(x, y, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, 1.0);
+      p4->Fill(x, y, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p3->Fill(x, y, 1.0);
+      p4->Fill(x, y, 1.0);
+   }
+
+   TList *list = new TList;
+   list->Add(p2);
+   list->Add(p3);
+
+   p1->Merge(list);
+
+   bool ret = equals("Merge1DProf", p1, p4, cmpOptStats, 1E-10);
+   delete p1;
+   delete p2;
+   delete p3;
+   return ret;
+}
+
+bool testMerge2D() 
+{
+   TH2D* h1 = new TH2D("merge2D-h1", "h1-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   TH2D* h2 = new TH2D("merge2D-h2", "h2-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   TH2D* h3 = new TH2D("merge2D-h3", "h3-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   TH2D* h4 = new TH2D("merge2D-h4", "h4-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+
+   h1->Sumw2();h2->Sumw2();h3->Sumw2();
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h1->Fill(x, y, 1.0);
+      h4->Fill(x, y, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h2->Fill(x, y, 1.0);
+      h4->Fill(x, y, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h3->Fill(x, y, 1.0);
+      h4->Fill(x, y, 1.0);
+   }
+
+   TList *list = new TList;
+   list->Add(h2);
+   list->Add(h3);
+
+   h1->Merge(list);
+
+   bool ret = equals("Merge2D", h1, h4, cmpOptStats, 1E-10);
+   delete h1;
+   delete h2;
+   delete h3;
+   return ret;
+}
+
+bool testMergeProf2D() 
+{
+   TProfile2D* p1 = new TProfile2D("merge2D-p1", "p1-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   TProfile2D* p2 = new TProfile2D("merge2D-p2", "p2-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   TProfile2D* p3 = new TProfile2D("merge2D-p3", "p3-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   TProfile2D* p4 = new TProfile2D("merge2D-p4", "p4-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, z, 1.0);
+      p4->Fill(x, y, z, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, z, 1.0);
+      p4->Fill(x, y, z, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p3->Fill(x, y, z, 1.0);
+      p4->Fill(x, y, z, 1.0);
+   }
+
+   TList *list = new TList;
+   list->Add(p2);
+   list->Add(p3);
+
+   p1->Merge(list);
+
+   bool ret = equals("Merge2DProf", p1, p4, cmpOptStats, 1E-10);
+   delete p1;
+   delete p2;
+   delete p3;
+   return ret;
 }
 
-// end stresHistRebin.cxx file
+bool testMerge3D() 
+{
+   TH3D* h1 = new TH3D("merge3D-h1", "h1-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   TH3D* h2 = new TH3D("merge3D-h2", "h2-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   TH3D* h3 = new TH3D("merge3D-h3", "h3-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
+   TH3D* h4 = new TH3D("merge3D-h4", "h4-Title",
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange,
+                       numberOfBins, minRange, maxRange);
 
-// old stresHistRebin.cxx file
+   h1->Sumw2();h2->Sumw2();h3->Sumw2();
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h1->Fill(x, y, z, 1.0);
+      h4->Fill(x, y, z, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h2->Fill(x, y, z, 1.0);
+      h4->Fill(x, y, z, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      h3->Fill(x, y, z, 1.0);
+      h4->Fill(x, y, z, 1.0);
+   }
+
+   TList *list = new TList;
+   list->Add(h2);
+   list->Add(h3);
+
+   h1->Merge(list);
+
+   bool ret = equals("Merge3D", h1, h4, cmpOptStats, 1E-10);
+   delete h1;
+   delete h2;
+   delete h3;
+   return ret;
+}
+
+bool testMergeProf3D() 
+{
+   TProfile3D* p1 = new TProfile3D("merge3D-p1", "p1-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   TProfile3D* p2 = new TProfile3D("merge3D-p2", "p2-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   TProfile3D* p3 = new TProfile3D("merge3D-p3", "p3-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+   TProfile3D* p4 = new TProfile3D("merge3D-p4", "p4-Title",
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange,
+                                   numberOfBins, minRange, maxRange);
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p1->Fill(x, y, z, t, 1.0);
+      p4->Fill(x, y, z, t, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p2->Fill(x, y, z, t, 1.0);
+      p4->Fill(x, y, z, t, 1.0);
+   }
+
+   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
+      Double_t x = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t y = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      Double_t t = r.Uniform(0.9 * minRange, 1.1 * maxRange);
+      p3->Fill(x, y, z, t, 1.0);
+      p4->Fill(x, y, z, t, 1.0);
+   }
+
+   TList *list = new TList;
+   list->Add(p2);
+   list->Add(p3);
+
+   p1->Merge(list);
+
+   bool ret = equals("Merge3DProf", p1, p4, cmpOptStats, 1E-10);
+   delete p1;
+   delete p2;
+   delete p3;
+   return ret;
+}
+
+bool testLabel()
+{
+   TH1D* h1 = new TH1D("lD1-h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("lD1-h2", "h2-Title", numberOfBins, minRange, maxRange);
+
+   for ( Int_t bin = 1; bin <= h1->GetNbinsX() ; ++bin ) {
+      ostringstream label;
+      label << bin;
+      h2->GetXaxis()->SetBinLabel(bin, label.str().c_str());
+   }
+
+   for ( Int_t e = 0; e < nEvents; ++e ) {
+      Double_t value = r.Uniform(minRange, maxRange);
+      Int_t bin = h1->GetXaxis()->FindBin(value);
+      h1->AddBinContent( bin, 1.0);
+
+      ostringstream label;
+      label << bin;
+      h2->Fill(label.str().c_str(), 1.0);
+   }
+
+   bool status = equals("Fill(char*)", h1, h2, cmpOptStats, 1E-13);
+   delete h1;
+   return status;
+}
 
 bool testIntegerRebin()
 {
@@ -1322,7 +1842,7 @@ bool testArrayRebin()
    std::sort(rebinArray, rebinArray + rebin);
    for ( Int_t i = 0; i < rebin; ++i ) {
       rebinArray[i] = TMath::Nint( rebinArray[i] * ( h1->GetNbinsX() - 2 ) + 2 );
-      rebinArray[i] = h1->GetBinLowEdge( h1->FindBin(rebinArray[i]) );
+      rebinArray[i] = h1->GetBinLowEdge( h1->GetXaxis()->FindBin( rebinArray[i] ) );
    }
    
 
@@ -1372,25 +1892,39 @@ bool test2DRebin()
    return equals("TestIntRebin2D", h2d2, h3, cmpOptStats);
 }
 
-bool stressHistRebin()
+bool testSparseRebin1() 
 {
-   const unsigned int numberOfTests = 4;
-   pointer2Test testPointer[numberOfTests] = { testIntegerRebin, 
-                                               testIntegerRebinNoName,
-                                               testArrayRebin,
-                                               test2DRebin };
-
-   bool status = false;
-   for ( unsigned int i = 0; i < numberOfTests; ++i )
-      status |= testPointer[i]();
-
-   return status;
+   const int rebin = 2;//TMath::Nint( r.Uniform(minRebin, maxRebin) );
+
+   Int_t bsizeRebin[] = { 2,//TMath::Nint( r.Uniform(1, 5) ),
+                          2,//TMath::Nint( r.Uniform(1, 5) ),
+                          2};//TMath::Nint( r.Uniform(1, 5) )};
+
+   Int_t bsize[] = { bsizeRebin[0] * rebin,
+                     bsizeRebin[1] * rebin,
+                     bsizeRebin[2] * rebin};
+                    
+   Double_t xmin[] = {minRange, minRange, minRange};
+   Double_t xmax[] = {maxRange, maxRange, maxRange};
+   THnSparseD* s1 = new THnSparseD("rebin1-s1","s1-Title", 3, bsize, xmin, xmax);
+   THnSparseD* s2 = new THnSparseD("rebin1-s2","s2-Title", 3, bsizeRebin, xmin, xmax);
+
+   for ( Int_t i = 0; i < nEvents; ++i ) {
+      Double_t points[3];
+      points[0] = r.Uniform( minRange * .9 , maxRange * 1.1 );
+      points[1] = r.Uniform( minRange * .9 , maxRange * 1.1 );
+      points[2] = r.Uniform( minRange * .9 , maxRange * 1.1 );
+      s1->Fill(points);
+      s2->Fill(points);
+   }
+
+   THnSparse* s3 = s1->Rebin(rebin);
+
+   bool ret = equals("THnSparse Rebin 1", s2, s3, cmpOptStats | cmpOptPrint );
+   delete s1;
+   return ret;
 }
 
-// end stressHistRebin.cxx file
- 
-
-// old stressHistProj file 
 
 // In case of deviation, the profiles' content will not work anymore
 // try only for testing the statistics
@@ -1400,7 +1934,9 @@ static const double centre_deviation = 0.0;
 class ProjectionTester {
    
 private:
-   static const unsigned int binsize = 10;
+   static const unsigned int binsizeX =  8;
+   static const unsigned int binsizeY = 10;
+   static const unsigned int binsizeZ = 12;
    static const int lower_limit = 0;
    static const int upper_limit = 10;
 
@@ -1463,59 +1999,59 @@ public:
    
    void CreateHistograms()
    {
-      h3 = new TH3D("h3","h3", binsize, lower_limit, upper_limit, 
-                               binsize, lower_limit, upper_limit, 
-                               binsize, lower_limit, upper_limit);
-
-      h2XY = new TH2D("h2XY", "h2XY", binsize, lower_limit, upper_limit, 
-                                      binsize, lower_limit, upper_limit);
-      h2XZ = new TH2D("h2XZ", "h2XZ", binsize, lower_limit, upper_limit, 
-                                      binsize, lower_limit, upper_limit);
-      h2YX = new TH2D("h2YX", "h2YX", binsize, lower_limit, upper_limit, 
-                                      binsize, lower_limit, upper_limit);
-      h2YZ = new TH2D("h2YZ", "h2YZ", binsize, lower_limit, upper_limit, 
-                                      binsize, lower_limit, upper_limit);
-      h2ZX = new TH2D("h2ZX", "h2ZX", binsize, lower_limit, upper_limit, 
-                                      binsize, lower_limit, upper_limit);
-      h2ZY = new TH2D("h2ZY", "h2ZY", binsize, lower_limit, upper_limit, 
-                                      binsize, lower_limit, upper_limit);
+      h3 = new TH3D("h3","h3", binsizeX, lower_limit, upper_limit, 
+                               binsizeY, lower_limit, upper_limit, 
+                               binsizeZ, lower_limit, upper_limit);
+
+      h2XY = new TH2D("h2XY", "h2XY", binsizeX, lower_limit, upper_limit, 
+                                      binsizeY, lower_limit, upper_limit);
+      h2XZ = new TH2D("h2XZ", "h2XZ", binsizeX, lower_limit, upper_limit, 
+                                      binsizeZ, lower_limit, upper_limit);
+      h2YX = new TH2D("h2YX", "h2YX", binsizeY, lower_limit, upper_limit, 
+                                      binsizeX, lower_limit, upper_limit);
+      h2YZ = new TH2D("h2YZ", "h2YZ", binsizeY, lower_limit, upper_limit, 
+                                      binsizeZ, lower_limit, upper_limit);
+      h2ZX = new TH2D("h2ZX", "h2ZX", binsizeZ, lower_limit, upper_limit, 
+                                      binsizeX, lower_limit, upper_limit);
+      h2ZY = new TH2D("h2ZY", "h2ZY", binsizeZ, lower_limit, upper_limit, 
+                                      binsizeY, lower_limit, upper_limit);
 
       // The bit is set for all the histograms (It's a statistic variable)
       TH1::StatOverflows(kTRUE);
 
-      h1X = new TH1D("h1X", "h1X", binsize, lower_limit, upper_limit);
-      h1Y = new TH1D("h1Y", "h1Y", binsize, lower_limit, upper_limit);
-      h1Z = new TH1D("h1Z", "h1Z", binsize, lower_limit, upper_limit);
-
-      h1XStats = new TH1D("h1XStats", "h1XStats", binsize, lower_limit, upper_limit);
-      h1YStats = new TH1D("h1YStats", "h1YStats", binsize, lower_limit, upper_limit);
-      h1ZStats = new TH1D("h1ZStats", "h1ZStats", binsize, lower_limit, upper_limit);
-
-      pe2XY = new TProfile2D("pe2XY", "pe2XY", binsize, lower_limit, upper_limit, 
-                                               binsize, lower_limit, upper_limit);
-      pe2XZ = new TProfile2D("pe2XZ", "pe2XZ", binsize, lower_limit, upper_limit, 
-                                               binsize, lower_limit, upper_limit);
-      pe2YX = new TProfile2D("pe2YX", "pe2YX", binsize, lower_limit, upper_limit, 
-                                               binsize, lower_limit, upper_limit);
-      pe2YZ = new TProfile2D("pe2YZ", "pe2YZ", binsize, lower_limit, upper_limit, 
-                                               binsize, lower_limit, upper_limit);
-      pe2ZX = new TProfile2D("pe2ZX", "pe2ZX", binsize, lower_limit, upper_limit, 
-                                               binsize, lower_limit, upper_limit);
-      pe2ZY = new TProfile2D("pe2ZY", "pe2ZY", binsize, lower_limit, upper_limit, 
-                                               binsize, lower_limit, upper_limit);
+      h1X = new TH1D("h1X", "h1X", binsizeX, lower_limit, upper_limit);
+      h1Y = new TH1D("h1Y", "h1Y", binsizeY, lower_limit, upper_limit);
+      h1Z = new TH1D("h1Z", "h1Z", binsizeZ, lower_limit, upper_limit);
+
+      h1XStats = new TH1D("h1XStats", "h1XStats", binsizeX, lower_limit, upper_limit);
+      h1YStats = new TH1D("h1YStats", "h1YStats", binsizeY, lower_limit, upper_limit);
+      h1ZStats = new TH1D("h1ZStats", "h1ZStats", binsizeZ, lower_limit, upper_limit);
+
+      pe2XY = new TProfile2D("pe2XY", "pe2XY", binsizeX, lower_limit, upper_limit, 
+                                               binsizeY, lower_limit, upper_limit);
+      pe2XZ = new TProfile2D("pe2XZ", "pe2XZ", binsizeX, lower_limit, upper_limit, 
+                                               binsizeZ, lower_limit, upper_limit);
+      pe2YX = new TProfile2D("pe2YX", "pe2YX", binsizeY, lower_limit, upper_limit, 
+                                               binsizeX, lower_limit, upper_limit);
+      pe2YZ = new TProfile2D("pe2YZ", "pe2YZ", binsizeY, lower_limit, upper_limit, 
+                                               binsizeZ, lower_limit, upper_limit);
+      pe2ZX = new TProfile2D("pe2ZX", "pe2ZX", binsizeZ, lower_limit, upper_limit, 
+                                               binsizeX, lower_limit, upper_limit);
+      pe2ZY = new TProfile2D("pe2ZY", "pe2ZY", binsizeZ, lower_limit, upper_limit, 
+                                               binsizeY, lower_limit, upper_limit);
       
-      h2wXY = new TH2D("h2wXY", "h2wXY", binsize, lower_limit, upper_limit, 
-                                         binsize, lower_limit, upper_limit);
-      h2wXZ = new TH2D("h2wXZ", "h2wXZ", binsize, lower_limit, upper_limit, 
-                                         binsize, lower_limit, upper_limit);
-      h2wYX = new TH2D("h2wYX", "h2wYX", binsize, lower_limit, upper_limit, 
-                                         binsize, lower_limit, upper_limit);
-      h2wYZ = new TH2D("h2wYZ", "h2wYZ", binsize, lower_limit, upper_limit, 
-                                         binsize, lower_limit, upper_limit);
-      h2wZX = new TH2D("h2wZX", "h2wZX", binsize, lower_limit, upper_limit, 
-                                         binsize, lower_limit, upper_limit);
-      h2wZY = new TH2D("h2wZY", "h2wZY", binsize, lower_limit, upper_limit, 
-                                         binsize, lower_limit, upper_limit);
+      h2wXY = new TH2D("h2wXY", "h2wXY", binsizeX, lower_limit, upper_limit, 
+                                         binsizeY, lower_limit, upper_limit);
+      h2wXZ = new TH2D("h2wXZ", "h2wXZ", binsizeX, lower_limit, upper_limit, 
+                                         binsizeZ, lower_limit, upper_limit);
+      h2wYX = new TH2D("h2wYX", "h2wYX", binsizeY, lower_limit, upper_limit, 
+                                         binsizeX, lower_limit, upper_limit);
+      h2wYZ = new TH2D("h2wYZ", "h2wYZ", binsizeY, lower_limit, upper_limit, 
+                                         binsizeZ, lower_limit, upper_limit);
+      h2wZX = new TH2D("h2wZX", "h2wZX", binsizeZ, lower_limit, upper_limit, 
+                                         binsizeX, lower_limit, upper_limit);
+      h2wZY = new TH2D("h2wZY", "h2wZY", binsizeZ, lower_limit, upper_limit, 
+                                         binsizeY, lower_limit, upper_limit);
 
       h2wXY->Sumw2();
       h2wXZ->Sumw2();
@@ -1524,19 +2060,19 @@ public:
       h2wZX->Sumw2();
       h2wZY->Sumw2();
 
-      pe1XY = new TProfile("pe1XY", "pe1XY", binsize, lower_limit, upper_limit);
-      pe1XZ = new TProfile("pe1XZ", "pe1XZ", binsize, lower_limit, upper_limit);
-      pe1YX = new TProfile("pe1YX", "pe1YX", binsize, lower_limit, upper_limit);
-      pe1YZ = new TProfile("pe1YZ", "pe1YZ", binsize, lower_limit, upper_limit);
-      pe1ZX = new TProfile("pe1ZX", "pe1ZX", binsize, lower_limit, upper_limit);
-      pe1ZY = new TProfile("pe1ZY", "pe1ZY", binsize, lower_limit, upper_limit);
+      pe1XY = new TProfile("pe1XY", "pe1XY", binsizeX, lower_limit, upper_limit);
+      pe1XZ = new TProfile("pe1XZ", "pe1XZ", binsizeX, lower_limit, upper_limit);
+      pe1YX = new TProfile("pe1YX", "pe1YX", binsizeY, lower_limit, upper_limit);
+      pe1YZ = new TProfile("pe1YZ", "pe1YZ", binsizeY, lower_limit, upper_limit);
+      pe1ZX = new TProfile("pe1ZX", "pe1ZX", binsizeZ, lower_limit, upper_limit);
+      pe1ZY = new TProfile("pe1ZY", "pe1ZY", binsizeZ, lower_limit, upper_limit);
 
-      hw1XY = new TH1D("hw1XY", "hw1XY", binsize, lower_limit, upper_limit);
-      hw1XZ = new TH1D("hw1XZ", "hw1XZ", binsize, lower_limit, upper_limit);
-      hw1YX = new TH1D("hw1YX", "hw1YX", binsize, lower_limit, upper_limit);
-      hw1YZ = new TH1D("hw1YZ", "hw1YZ", binsize, lower_limit, upper_limit);
-      hw1ZX = new TH1D("hw1ZX", "hw1ZX", binsize, lower_limit, upper_limit);
-      hw1ZY = new TH1D("hw1ZY", "hw1ZY", binsize, lower_limit, upper_limit);
+      hw1XY = new TH1D("hw1XY", "hw1XY", binsizeX, lower_limit, upper_limit);
+      hw1XZ = new TH1D("hw1XZ", "hw1XZ", binsizeX, lower_limit, upper_limit);
+      hw1YX = new TH1D("hw1YX", "hw1YX", binsizeY, lower_limit, upper_limit);
+      hw1YZ = new TH1D("hw1YZ", "hw1YZ", binsizeY, lower_limit, upper_limit);
+      hw1ZX = new TH1D("hw1ZX", "hw1ZX", binsizeZ, lower_limit, upper_limit);
+      hw1ZY = new TH1D("hw1ZY", "hw1ZY", binsizeZ, lower_limit, upper_limit);
 
       hw1XZ->Sumw2();
       hw1XY->Sumw2();
@@ -1545,7 +2081,7 @@ public:
       hw1ZX->Sumw2();
       hw1ZY->Sumw2();
 
-      Int_t bsize[] = {binsize, binsize, binsize};
+      Int_t bsize[] = {binsizeX, binsizeY, binsizeZ};
       Double_t xmin[] = {lower_limit, lower_limit, lower_limit};
       Double_t xmax[] = {upper_limit, upper_limit, upper_limit};
       s3 = new THnSparseD("s3","s3", 3, bsize, xmin, xmax);
@@ -1768,7 +2304,8 @@ public:
                        h3->GetYaxis()->FindBin(y) >= ymin && h3->GetYaxis()->FindBin(y) <= ymax &&
                        h3->GetZaxis()->FindBin(z) >= zmin && h3->GetZaxis()->FindBin(z) <= zmax )
                   {
-                     cout << "Filling (" << x << "," << y << "," << z << ")!" << endl;
+                     if ( defaultEqualOptions & cmpOptPrint )
+                        cout << "Filling (" << x << "," << y << "," << z << ")!" << endl;
                      
                      h2XY->Fill(x,y);
                      h2XZ->Fill(x,z);
@@ -1863,7 +2400,8 @@ public:
       status += equals("TH3 -> ZX", h2ZX, (TH2D*) h3->Project3D("XZ"), options);
       status += equals("TH3 -> ZY", h2ZY, (TH2D*) h3->Project3D("YZ"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // TH1 derived from TH3
       options = cmpOptStats;
@@ -1871,7 +2409,8 @@ public:
       status += equals("TH3 -> Y", h1Y, (TH1D*) h3->Project3D("y"), options);
       status += equals("TH3 -> Z", h1Z, (TH1D*) h3->Project3D("z"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // TH1 derived from h2XY
       options = cmpOptStats;
@@ -1893,7 +2432,8 @@ public:
       status += equals("TH2ZY -> Z", h1Z, (TH1D*) h2ZY->ProjectionX("z"), options);
       status += equals("TH2ZY -> Y", h1Y, (TH1D*) h2ZY->ProjectionY("y"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // Now the histograms comming from the Profiles!
       options = cmpOptStats;
@@ -1904,7 +2444,8 @@ public:
       status += equals("TH3 -> PBZX", h2ZX, (TH2D*) h3->Project3DProfile("xz UF OF")->ProjectionXY("5", "B"), options);
       status += equals("TH3 -> PBZY", h2ZY, (TH2D*) h3->Project3DProfile("yz UF OF")->ProjectionXY("6", "B"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // test directly project3dprofile
       options = cmpOptStats;
@@ -1915,7 +2456,8 @@ public:
       status += equals("TH3 -> PZX", (TH2D*) pe2ZX, (TH2D*) h3->Project3DProfile("xz  UF OF"), options);
       status += equals("TH3 -> PZY", (TH2D*) pe2ZY, (TH2D*) h3->Project3DProfile("yz  UF OF"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // test option E of ProjectionXY
       options = 0;
@@ -1926,7 +2468,8 @@ public:
       status += equals("TH3 -> PEZX", (TH2D*) pe2ZX, (TH2D*) h3->Project3DProfile("xz  UF OF")->ProjectionXY("5", "E"), options);
       status += equals("TH3 -> PEZY", (TH2D*) pe2ZY, (TH2D*) h3->Project3DProfile("yz  UF OF")->ProjectionXY("6", "E"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // test option W of ProjectionXY
       
@@ -1939,100 +2482,109 @@ public:
       status += equals("TH3 -> PWZX", (TH2D*) h2wZX, (TH2D*) h3->Project3DProfile("xz  UF OF")->ProjectionXY("5", "W"), options);
       status += equals("TH3 -> PWZY", (TH2D*) h2wZY, (TH2D*) h3->Project3DProfile("yz  UF OF")->ProjectionXY("6", "W"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // test 1D histograms
       options = cmpOptStats;
-      status += equals("TH2XY -> PBX", h1X, (TH1D*) h2XY->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2XY -> PBX", h1Y, (TH1D*) h2XY->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2XZ -> PBX", h1X, (TH1D*) h2XZ->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2XZ -> PBZ", h1Z, (TH1D*) h2XZ->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2YX -> PBY", h1Y, (TH1D*) h2YX->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2YX -> PBX", h1X, (TH1D*) h2YX->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2YZ -> PBY", h1Y, (TH1D*) h2YZ->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2YZ -> PBZ", h1Z, (TH1D*) h2YZ->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2ZX -> PBZ", h1Z, (TH1D*) h2ZX->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2ZX -> PBX", h1X, (TH1D*) h2ZX->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2ZY -> PBZ", h1Z, (TH1D*) h2ZY->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
-      status += equals("TH2ZY -> PBY", h1Y, (TH1D*) h2ZY->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      // ProfileX re-use the same histo if sme name is given. 
+      // need to give a diffrent name for each projectino (x,y,Z) otherwise we end-up in different bins
+      // t.b.d: ProfileX make a new histo if non compatible
+      status += equals("TH2XY -> PBX", h1X, (TH1D*) h2XY->ProfileX("PBX", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      status += equals("TH2XY -> PBY", h1Y, (TH1D*) h2XY->ProfileY("PBY", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options );
+      status += equals("TH2XZ -> PBX", h1X, (TH1D*) h2XZ->ProfileX("PBX", 0,h2XZ->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      status += equals("TH2XZ -> PBZ", h1Z, (TH1D*) h2XZ->ProfileY("PBZ", 0,h2XZ->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options,1E-12);
+      status += equals("TH2YX -> PBY", h1Y, (TH1D*) h2YX->ProfileX("PBY", 0,h2YX->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      status += equals("TH2YX -> PBX", h1X, (TH1D*) h2YX->ProfileY("PBX", 0,h2YX->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      status += equals("TH2YZ -> PBY", h1Y, (TH1D*) h2YZ->ProfileX("PBY", 0,h2YZ->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      status += equals("TH2YZ -> PBZ", h1Z, (TH1D*) h2YZ->ProfileY("PBZ", 0,h2YZ->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options,1E-12);
+      status += equals("TH2ZX -> PBZ", h1Z, (TH1D*) h2ZX->ProfileX("PBZ", 0,h2ZX->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options,1E-12);
+      status += equals("TH2ZX -> PBX", h1X, (TH1D*) h2ZX->ProfileY("PBX", 0,h2ZX->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
+      status += equals("TH2ZY -> PBZ", h1Z, (TH1D*) h2ZY->ProfileX("PBZ", 0,h2ZY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "B"),options,1E-12);
+      status += equals("TH2ZY -> PBY", h1Y, (TH1D*) h2ZY->ProfileY("PBY", 0,h2ZY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "B"),options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
 
       // 1D testing direct profiles 
       options = cmpOptStats;
-      status += equals("TH2XY -> PX", pe1XY, (TH1D*) h2XY->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1), options);
-      status += equals("TH2XY -> PX", pe1YX, (TH1D*) h2XY->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1), options);
-      status += equals("TH2XZ -> PX", pe1XZ, (TH1D*) h2XZ->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1), options);
-      status += equals("TH2XZ -> PZ", pe1ZX, (TH1D*) h2XZ->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1), options);
-      status += equals("TH2YX -> PY", pe1YX, (TH1D*) h2YX->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1), options);
-      status += equals("TH2YX -> PX", pe1XY, (TH1D*) h2YX->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1), options);
-      status += equals("TH2YZ -> PY", pe1YZ, (TH1D*) h2YZ->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1), options);
-      status += equals("TH2YZ -> PZ", pe1ZY, (TH1D*) h2YZ->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1), options);
-      status += equals("TH2ZX -> PZ", pe1ZX, (TH1D*) h2ZX->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1), options);
-      status += equals("TH2ZX -> PX", pe1XZ, (TH1D*) h2ZX->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1), options);
-      status += equals("TH2ZY -> PZ", pe1ZY, (TH1D*) h2ZY->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1), options);
-      status += equals("TH2ZY -> PY", pe1YZ, (TH1D*) h2ZY->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2XY -> PX", pe1XY, (TH1D*) h2XY->ProfileX("PX", 0,h2XY->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2XY -> PY", pe1YX, (TH1D*) h2XY->ProfileY("PY", 0,h2XY->GetXaxis()->GetNbins()+1), options);
+      status += equals("TH2XZ -> PX", pe1XZ, (TH1D*) h2XZ->ProfileX("PX", 0,h2XZ->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2XZ -> PZ", pe1ZX, (TH1D*) h2XZ->ProfileY("PZ", 0,h2XZ->GetXaxis()->GetNbins()+1), options);
+      status += equals("TH2YX -> PY", pe1YX, (TH1D*) h2YX->ProfileX("PY", 0,h2YX->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2YX -> PX", pe1XY, (TH1D*) h2YX->ProfileY("PX", 0,h2YX->GetXaxis()->GetNbins()+1), options);
+      status += equals("TH2YZ -> PY", pe1YZ, (TH1D*) h2YZ->ProfileX("PY", 0,h2YZ->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2YZ -> PZ", pe1ZY, (TH1D*) h2YZ->ProfileY("PZ", 0,h2YZ->GetXaxis()->GetNbins()+1), options);
+      status += equals("TH2ZX -> PZ", pe1ZX, (TH1D*) h2ZX->ProfileX("PZ", 0,h2ZX->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2ZX -> PX", pe1XZ, (TH1D*) h2ZX->ProfileY("PX", 0,h2ZX->GetXaxis()->GetNbins()+1), options);
+      status += equals("TH2ZY -> PZ", pe1ZY, (TH1D*) h2ZY->ProfileX("PZ", 0,h2ZY->GetYaxis()->GetNbins()+1), options);
+      status += equals("TH2ZY -> PY", pe1YZ, (TH1D*) h2ZY->ProfileY("PY", 0,h2ZY->GetXaxis()->GetNbins()+1), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
 
       // 1D testing e profiles
       options = 0;
       status += equals("TH2XY -> PEX", pe1XY, 
-                       (TH1D*) h2XY->ProfileX("8", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
-      status += equals("TH2XY -> PEX", pe1YX, 
-                       (TH1D*) h2XY->ProfileY("8", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2XY->ProfileX("PEX", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+      status += equals("TH2XY -> PEY", pe1YX, 
+                       (TH1D*) h2XY->ProfileY("PEY", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2XZ -> PEX", pe1XZ, 
-                       (TH1D*) h2XZ->ProfileX("8", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2XZ->ProfileX("PEX", 0,h2XZ->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2XZ -> PEZ", pe1ZX, 
-                       (TH1D*) h2XZ->ProfileY("8", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2XZ->ProfileY("PEZ", 0,h2XZ->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2YX -> PEY", pe1YX, 
-                       (TH1D*) h2YX->ProfileX("8", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2YX->ProfileX("PEY", 0,h2YX->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2YX -> PEX", pe1XY, 
-                       (TH1D*) h2YX->ProfileY("8", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2YX->ProfileY("PEX", 0,h2YX->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2YZ -> PEY", pe1YZ, 
-                       (TH1D*) h2YZ->ProfileX("8", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2YZ->ProfileX("PEY", 0,h2YZ->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2YZ -> PEZ", pe1ZY, 
-                       (TH1D*) h2YZ->ProfileY("8", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2YZ->ProfileY("PEZ", 0,h2YZ->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2ZX -> PEZ", pe1ZX, 
-                       (TH1D*) h2ZX->ProfileX("8", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2ZX->ProfileX("PEZ", 0,h2ZX->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2ZX -> PEX", pe1XZ, 
-                       (TH1D*) h2ZX->ProfileY("8", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2ZX->ProfileY("PEX", 0,h2ZX->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2ZY -> PEZ", pe1ZY, 
-                       (TH1D*) h2ZY->ProfileX("8", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2ZY->ProfileX("PEZ", 0,h2ZY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       status += equals("TH2ZY -> PEY", pe1YZ, 
-                       (TH1D*) h2ZY->ProfileY("8", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
+                       (TH1D*) h2ZY->ProfileY("PEY", 0,h2ZY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "E"), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
 
       // 1D testing w profiles
       // The error is not properly propagated when build with weights :S
       if ( buildWithWeights ) options = cmpOptNoError;
       status += equals("TH2XY -> PWX", hw1XY, 
-                       (TH1D*) h2XY->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
-      status += equals("TH2XY -> PWX", hw1YX, 
-                       (TH1D*) h2XY->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2XY->ProfileX("PWX", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+      status += equals("TH2XY -> PWY", hw1YX, 
+                       (TH1D*) h2XY->ProfileY("PWY", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2XZ -> PWX", hw1XZ, 
-                       (TH1D*) h2XZ->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2XZ->ProfileX("PWX", 0,h2XZ->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2XZ -> PWZ", hw1ZX, 
-                       (TH1D*) h2XZ->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2XZ->ProfileY("PWZ", 0,h2XZ->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2YX -> PWY", hw1YX, 
-                       (TH1D*) h2YX->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2YX->ProfileX("PWY", 0,h2YX->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2YX -> PWX", hw1XY, 
-                       (TH1D*) h2YX->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2YX->ProfileY("PWX", 0,h2YX->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2YZ -> PWY", hw1YZ, 
-                       (TH1D*) h2YZ->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2YZ->ProfileX("PWY", 0,h2YZ->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2YZ -> PWZ", hw1ZY, 
-                       (TH1D*) h2YZ->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2YZ->ProfileY("PWZ", 0,h2YZ->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2ZX -> PWZ", hw1ZX, 
-                       (TH1D*) h2ZX->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2ZX->ProfileX("PWZ", 0,h2ZX->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2ZX -> PWX", hw1XZ, 
-                       (TH1D*) h2ZX->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2ZX->ProfileY("PWX", 0,h2ZX->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2ZY -> PWZ", hw1ZY, 
-                       (TH1D*) h2ZY->ProfileX("7", 0,h2XY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2ZY->ProfileX("PWZ", 0,h2ZY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
       status += equals("TH2ZY -> PWY", hw1YZ, 
-                       (TH1D*) h2ZY->ProfileY("7", 0,h2XY->GetYaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+                       (TH1D*) h2ZY->ProfileY("PWY", 0,h2ZY->GetXaxis()->GetNbins()+1)->ProjectionX("1", "W"), options);
+
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
       
       // TH2 derived from STH3
       options = cmpOptStats;
@@ -2040,10 +2592,11 @@ public:
       status += equals("STH3 -> XZ", h2XZ, (TH2D*) s3->Projection(2,0), options);
       status += equals("STH3 -> YX", h2YX, (TH2D*) s3->Projection(0,1), options);
       status += equals("STH3 -> YZ", h2YZ, (TH2D*) s3->Projection(2,1), options);
-      status += equals("STH3 -> ZX", h2ZX, (TH2D*) s3->Projection(0,2), options);
-      status += equals("STH3 -> ZY", h2ZY, (TH2D*) s3->Projection(1,2), options);
+      status += equals("STH3 -> ZX", h2ZX, (TH2D*) s3->Projection(0,2), options); 
+      status += equals("STH3 -> ZY", h2ZY, (TH2D*) s3->Projection(1,2), options); 
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
 
       // TH1 derived from STH3
       options = cmpOptStats;
@@ -2051,49 +2604,14 @@ public:
       status += equals("STH3 -> Y", h1Y, (TH1D*) s3->Projection(1), options);
       status += equals("STH3 -> Z", h1Z, (TH1D*) s3->Projection(2), options);
       options = 0;
-      cout << "----------------------------------------------" << endl;
+      if ( defaultEqualOptions & cmpOptPrint )
+         cout << "----------------------------------------------" << endl;
 
       return status;
    }
    
 };
 
-int stressHistProj(bool testWithoutWeights = true,
-                   bool testWithWeights = true)
-{
-   int status = 0;
-   
-   if ( testWithoutWeights )
-   {
-      cout << "**********************************\n"
-           << "       Test without weights       \n" 
-           << "**********************************\n"
-           << endl;
-      
-      ProjectionTester ht;
-      ht.buildHistograms();
-      //ht.buildHistograms(2,4,5,6,8,10);
-      status += ht.compareHistograms();
-   }
-
-   if ( testWithWeights )
-   {
-      cout << "**********************************\n"
-           << "        Test with weights         \n" 
-           << "**********************************\n"
-           << endl;
-
-      ProjectionTester ht;
-      ht.buildHistogramsWithWeights();
-      status += ht.compareHistograms();
-   }
-
-   return status;
-}
-
-
-// end of stressHistProj file
-
 int main(int argc, char** argv)
 {
    r.SetSeed(0);
@@ -2103,31 +2621,161 @@ int main(int argc, char** argv)
    if ( __DRAW__ )
       theApp = new TApplication("App",&argc,argv);
 
-   bool GlobalStatus = false;
-   bool status = false;
+   int GlobalStatus = false;
+   int status = false;
 
-   ostringstream output;
-   output << "\nTEST RESULTS\n\n";
+   cout << "****************************************************************************" <<endl;
+   cout << "*  Starting  stress  H I S T O G R A M                                     *" <<endl;
+   cout << "****************************************************************************" <<endl;
 
-   cout << "\nstressHistProj\n" << endl;
-   status = stressHistProj();
-   GlobalStatus |= status;
-   output << "stressHistProj Test.............." 
-          << (status?"FAILED":"OK") << endl;
+   // Test 1
+   if ( defaultEqualOptions & cmpOptPrint )
+      cout << "**********************************\n"
+           << "       Test without weights       \n" 
+           << "**********************************\n"
+           << endl;
+   
+   ProjectionTester* ht = new ProjectionTester();
+   ht->buildHistograms();
+   //ht->buildHistograms(2,4,5,6,8,10);
+   status = ht->compareHistograms();
+   GlobalStatus += status;
+   printResult("Testing Projections without weights..............................", status);
+   delete ht;
+
+   // Test 2
+   if ( defaultEqualOptions & cmpOptPrint )
+      cout << "**********************************\n"
+           << "        Test with weights         \n" 
+           << "**********************************\n"
+           << endl;
+   
+   ProjectionTester* ht2 = new ProjectionTester();
+   ht2->buildHistogramsWithWeights();
+   status = ht2->compareHistograms();
+   GlobalStatus += status;
+   printResult("Testing Projections with weights.................................", status);
+   delete ht2;
+   
+   // Test 3
+   const unsigned int numberOfRebin = 4;
+   pointer2Test rebinTestPointer[numberOfRebin] = { testIntegerRebin, 
+                                                    testIntegerRebinNoName,
+                                                    testArrayRebin,
+                                                    test2DRebin,
+                                                    /*testSparseRebin1*/};
+   struct TTestSuite rebinTestSuite = { numberOfRebin, 
+                                        "Histogram Rebinning..............................................",
+                                        rebinTestPointer };
+
+   // testSparseRebin1 fails. To be looked. Also to be done the
+   // testSparseRebin2, but this is calling the first method, so also
+   // frozen until fixed.
+
+   // Test 4
+   // Add Tests
+   const unsigned int numberOfAdds = 12;
+   pointer2Test addTestPointer[numberOfAdds] = { testAdd1,    testAddProfile1, 
+                                                 testAdd2,    testAddProfile2,
+                                                 testAdd2D1,  testAdd2DProfile1,
+                                                 testAdd2D2,  testAdd2DProfile2,
+                                                 testAdd3D1,  testAdd3DProfile1,
+                                                 testAdd3D2,  testAdd3DProfile2
+   };
+   struct TTestSuite addTestSuite = { numberOfAdds, 
+                                      "Add tests for 1D, 2D and 3D Histograms and Profiles..............",
+                                      addTestPointer };
+
+   // Test 5
+   // Multiply Tests
+   const unsigned int numberOfMultiply = 6;
+   pointer2Test multiplyTestPointer[numberOfMultiply] = { testMul1,    testMul2,
+                                                          testMul2D1,  testMul2D2,
+                                                          testMul3D1,  testMul3D2
+   };
+   struct TTestSuite multiplyTestSuite = { numberOfMultiply, 
+                                           "Multiply tests for 1D, 2D and 3D Histograms......................",
+                                           multiplyTestPointer };
 
-   cout << "\nstressHistRebin\n" << endl;
-   status = stressHistRebin();
-   GlobalStatus |= status;
-   output << "stressHistRebin Test............."
-          << (status?"FAILED":"OK") << endl;
+   // Still to do: testDivide2, testDivide2D1, testDivide2D2 and
+   // testDivide3D1, testDivide3D2. 
 
-   cout << "\nstressHistOpts\n" << endl;
-   status = stressHistOpts();
-   GlobalStatus |= status;
-   output << "stressHistOpts Test.............."
-          << (status?"FAILED":"OK") << endl;
+   // It depends on whether we can solve the problem with the 1D test
+   // already done. It seems like there is something wrong with the
+   // way the errors are being calculated. We have a formula to
+   // calculate it by hand. Nevertheless the Divide method errors and
+   // the ones calculated by ourselves differ a bit from those (in the
+   // order of 1E-1).
 
-   cout << output.str() << endl;
+   // Test 6
+   // Copy Tests
+   const unsigned int numberOfCopy = 18;
+   pointer2Test copyTestPointer[numberOfCopy] = { testAssign1D,          testAssignProfile1D, 
+                                                  testCopyConstructor1D, testCopyConstructorProfile1D, 
+                                                  testClone1D,           testCloneProfile1D,
+                                                  testAssign2D,          testAssignProfile2D,
+                                                  testCopyConstructor2D, testCopyConstructorProfile2D,
+                                                  testClone2D,           testCloneProfile2D,
+                                                  testAssign3D,          testAssignProfile3D,
+                                                  testCopyConstructor3D, testCopyConstructorProfile3D,
+                                                  testClone3D,           testCloneProfile3D,
+   };
+   struct TTestSuite copyTestSuite = { numberOfCopy, 
+                                       "Copy tests for 1D, 2D and 3D Histograms and Profiles.............",
+                                       copyTestPointer };
+
+   // Test 7
+   // Readwrite Tests
+   const unsigned int numberOfReadwrite = 6;
+   pointer2Test readwriteTestPointer[numberOfReadwrite] = { testWriteRead1D,  testWriteReadProfile1D,
+                                                            testWriteRead2D,  testWriteReadProfile2D,
+                                                            testWriteRead3D,  testWriteReadProfile3D, 
+   };
+   struct TTestSuite readwriteTestSuite = { numberOfReadwrite, 
+                                            "Read/Write tests for 1D, 2D and 3D Histograms and Profiles.......",
+                                            readwriteTestPointer };
+
+   // Test 8
+   // Merge Tests
+   const unsigned int numberOfMerge = 6;
+   pointer2Test mergeTestPointer[numberOfMerge] = { testMerge1D,  testMergeProf1D,
+                                                    testMerge2D,  testMergeProf2D,
+                                                    testMerge3D,  testMergeProf3D
+   };
+   struct TTestSuite mergeTestSuite = { numberOfMerge, 
+                                        "Merge tests for 1D, 2D and 3D Histograms and Profiles............",
+                                        mergeTestPointer };
+   // Test 9
+   // Label Tests
+   const unsigned int numberOfLabel = 1;
+   pointer2Test labelTestPointer[numberOfLabel] = { testLabel
+   };
+   struct TTestSuite labelTestSuite = { numberOfLabel, 
+                                        "Label tests for 1D Histograms (TAxis)............................",
+                                        labelTestPointer };
+
+
+   // Combination of tests
+   const unsigned int numberOfSuits = 7;
+   struct TTestSuite* testSuite[numberOfSuits];
+   testSuite[0] = &rebinTestSuite;
+   testSuite[1] = &addTestSuite;
+   testSuite[2] = &multiplyTestSuite;
+   testSuite[3] = &copyTestSuite;
+   testSuite[4] = &readwriteTestSuite;
+   testSuite[5] = &mergeTestSuite;
+   testSuite[6] = &labelTestSuite;
+
+   status = 0;
+   for ( unsigned int i = 0; i < numberOfSuits; ++i ) {
+      bool internalStatus = false;
+      for ( unsigned int j = 0; j < testSuite[i]->nTests; ++j ) {
+         internalStatus |= testSuite[i]->tests[j]();
+      }
+      printResult( testSuite[i]->suiteName, internalStatus);
+      status += internalStatus;
+   }
+   GlobalStatus += status;
 
    if ( __DRAW__ ) {
       theApp->Run();
@@ -2135,7 +2783,7 @@ int main(int argc, char** argv)
       theApp = 0;
    }
 
-   return status;
+   return GlobalStatus;
 }
 
 ostream& operator<<(ostream& out, TH1D* h)
@@ -2148,15 +2796,78 @@ ostream& operator<<(ostream& out, TH1D* h)
    return out;
 }
 
+void printResult(const char* msg, bool status)
+{
+   static int counter = 1;
+   cout << "Test ";
+   cout.width(2);
+   cout<< counter << ": "
+       << msg
+       << (status?"FAILED":"OK") << endl;
+   counter += 1;
+}
+
 // Methods for histogram comparisions
 
+int equals(const char* msg, THnSparse* h1, THnSparse* h2, int options, double ERRORLIMIT)
+{
+   bool debug = options & cmpOptDebug;
+   bool compareError = ! (options & cmpOptNoError);
+   
+   int differents = 0;
+   
+   for ( int i = 0; i <= h1->GetAxis(0)->GetNbins() + 1; ++i )
+      for ( int j = 0; j <= h1->GetAxis(1)->GetNbins() + 1; ++j )
+         for ( int h = 0; h <= h1->GetAxis(2)->GetNbins() + 1; ++h )
+         {
+            Double_t x = h1->GetAxis(0)->GetBinCenter(i);
+            Double_t y = h1->GetAxis(1)->GetBinCenter(j);
+            Double_t z = h1->GetAxis(2)->GetBinCenter(h);
+            
+            Int_t bin[3] = {i, j, h};
+            
+            if (debug) {
+               cout << equals(x, h2->GetAxis(0)->GetBinCenter(i), ERRORLIMIT) << " "
+                    << equals(y, h2->GetAxis(1)->GetBinCenter(j), ERRORLIMIT) << " "
+                    << equals(z, h2->GetAxis(2)->GetBinCenter(h), ERRORLIMIT) << " "
+                    << "[" << x << "," << y << "," << z << "]: " 
+                    << h1->GetBinContent(bin) << " +/- " << h1->GetBinError(bin) << " | "
+                    << h2->GetBinContent(bin) << " +/- " << h2->GetBinError(bin)
+                    << " | " << equals(h1->GetBinContent(bin), h2->GetBinContent(bin), ERRORLIMIT)
+                    << " "   << equals(h1->GetBinError(bin)  , h2->GetBinError(bin),   ERRORLIMIT)
+                    << " "   << differents
+                    << " "   << (fabs(h1->GetBinContent(bin) - h2->GetBinContent(bin)))
+                    << endl;
+            }
+            differents += equals(x, h2->GetAxis(0)->GetBinCenter(i), ERRORLIMIT);
+            differents += equals(y, h2->GetAxis(1)->GetBinCenter(j), ERRORLIMIT);
+            differents += equals(z, h2->GetAxis(2)->GetBinCenter(h), ERRORLIMIT);
+            differents += equals(h1->GetBinContent(bin), h2->GetBinContent(bin), ERRORLIMIT);
+            if ( compareError )
+               differents += equals(h1->GetBinError(bin)  , h2->GetBinError(bin), ERRORLIMIT);
+         }
+   
+   // Statistical tests:
+   // No statistical tests possible for THnSparse so far...
+//    if ( compareStats )
+//       differents += compareStatistics( h1, h2, debug, ERRORLIMIT);
+   
+   cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
+   
+   delete h2;
+   
+   return differents;
+}
+
 int equals(const char* msg, TH3D* h1, TH3D* h2, int options, double ERRORLIMIT)
 {
+   options = options | defaultEqualOptions;
+   bool print = options & cmpOptPrint;
    bool debug = options & cmpOptDebug;
    bool compareError = ! (options & cmpOptNoError);
    bool compareStats = options & cmpOptStats;
    
-   bool differents = ( h1 == h2 ); // Check they are not the same histogram!
+   int differents = ( h1 == h2 ); // Check they are not the same histogram!
    if (debug) {
       cout << static_cast<void*>(h1) << " " << static_cast<void*>(h2) << " "
            << (h1 == h2 ) << " " << differents << endl;
@@ -2184,19 +2895,19 @@ int equals(const char* msg, TH3D* h1, TH3D* h2, int options, double ERRORLIMIT)
                  << " "   << (fabs(h1->GetBinContent(i,j,h) - h2->GetBinContent(i,j,h)))
                  << endl;
          }
-         differents |= (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
-         differents |= (bool) equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT);
-         differents |= (bool) equals(z, h2->GetZaxis()->GetBinCenter(h), ERRORLIMIT);
-         differents |= (bool) equals(h1->GetBinContent(i,j,h), h2->GetBinContent(i,j,h), ERRORLIMIT);
+         differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
+         differents += (bool) equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT);
+         differents += (bool) equals(z, h2->GetZaxis()->GetBinCenter(h), ERRORLIMIT);
+         differents += (bool) equals(h1->GetBinContent(i,j,h), h2->GetBinContent(i,j,h), ERRORLIMIT);
          if ( compareError )
-            differents |= (bool) equals(h1->GetBinError(i,j,h)  , h2->GetBinError(i,j,h), ERRORLIMIT);
+            differents += (bool) equals(h1->GetBinError(i,j,h)  , h2->GetBinError(i,j,h), ERRORLIMIT);
       }
    
    // Statistical tests:
    if ( compareStats )
-     differents |= (bool) compareStatistics( h1, h2, debug, ERRORLIMIT);
+      differents += compareStatistics( h1, h2, debug, ERRORLIMIT);
    
-   cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
+   if ( print || debug ) cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
    
    delete h2;
    
@@ -2205,11 +2916,13 @@ int equals(const char* msg, TH3D* h1, TH3D* h2, int options, double ERRORLIMIT)
 
 int equals(const char* msg, TH2D* h1, TH2D* h2, int options, double ERRORLIMIT)
 {
+   options = options | defaultEqualOptions;
+   bool print = options & cmpOptPrint;
    bool debug = options & cmpOptDebug;
    bool compareError = ! (options & cmpOptNoError);
    bool compareStats = options & cmpOptStats;
    
-   bool differents = ( h1 == h2 ); // Check they are not the same histogram!
+   int differents = ( h1 == h2 ); // Check they are not the same histogram!
    if (debug) {
       cout << static_cast<void*>(h1) << " " << static_cast<void*>(h2) << " "
            << (h1 == h2 ) << " " << differents << endl;
@@ -2234,18 +2947,18 @@ int equals(const char* msg, TH2D* h1, TH2D* h2, int options, double ERRORLIMIT)
                  << " "   << (fabs(h1->GetBinContent(i,j) - h2->GetBinContent(i,j)))
                  << endl;
          }
-         differents |= (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
-         differents |= (bool) equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT);
-         differents |= (bool) equals(h1->GetBinContent(i,j), h2->GetBinContent(i,j), ERRORLIMIT);
+         differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
+         differents += (bool) equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT);
+         differents += (bool) equals(h1->GetBinContent(i,j), h2->GetBinContent(i,j), ERRORLIMIT);
          if ( compareError )
-            differents |= (bool) equals(h1->GetBinError(i,j)  , h2->GetBinError(i,j), ERRORLIMIT);
+            differents += (bool) equals(h1->GetBinError(i,j)  , h2->GetBinError(i,j), ERRORLIMIT);
       }
    
    // Statistical tests:
    if ( compareStats )
-     differents |= (bool) compareStatistics( h1, h2, debug, ERRORLIMIT);
+      differents += compareStatistics( h1, h2, debug, ERRORLIMIT);
    
-   cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
+   if ( print || debug ) cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
    
    delete h2;
    
@@ -2254,11 +2967,18 @@ int equals(const char* msg, TH2D* h1, TH2D* h2, int options, double ERRORLIMIT)
 
 int equals(const char* msg, TH1D* h1, TH1D* h2, int options, double ERRORLIMIT)
 {
+   options = options | defaultEqualOptions;
+   bool print = options & cmpOptPrint;
    bool debug = options & cmpOptDebug;
    bool compareError = ! (options & cmpOptNoError);
    bool compareStats = options & cmpOptStats;
    
-   bool differents = ( h1 == h2 ); // Check they are not the same histogram!
+
+   if (debug) { 
+      cout << "Nbins  = " << h1->GetXaxis()->GetNbins() << " ,  " <<  h2->GetXaxis()->GetNbins() << endl;
+   }
+
+   int differents = ( h1 == h2 ); // Check they are not the same histogram!
    if (debug) {
       cout << static_cast<void*>(h1) << " " << static_cast<void*>(h2) << " "
            << (h1 == h2 ) << " " << differents << endl;
@@ -2278,18 +2998,18 @@ int equals(const char* msg, TH1D* h1, TH1D* h2, int options, double ERRORLIMIT)
               << " "   << differents
               << endl;
       }
-      differents |= (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
-      differents |= (bool) equals(h1->GetBinContent(i), h2->GetBinContent(i), ERRORLIMIT);
+      differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
+      differents += (bool) equals(h1->GetBinContent(i), h2->GetBinContent(i), ERRORLIMIT);
       
       if ( compareError )
-         differents |= (bool) equals(h1->GetBinError(i),   h2->GetBinError(i), ERRORLIMIT);
+         differents += (bool) equals(h1->GetBinError(i),   h2->GetBinError(i), ERRORLIMIT);
    }
    
    // Statistical tests:
    if ( compareStats )
-     differents |= (bool) compareStatistics( h1, h2, debug, ERRORLIMIT);
+      differents += compareStatistics( h1, h2, debug, ERRORLIMIT);
    
-   cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
+   if ( print || debug ) cout << msg << ": \t" << (differents?"FAILED":"OK") << endl;
    
    delete h2;
    
@@ -2303,7 +3023,7 @@ int equals(Double_t n1, Double_t n2, double ERRORLIMIT)
 
 int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT)
 {
-   bool differents = 0;
+   int differents = 0;
 
    int precLevel = gErrorIgnoreLevel; 
    // switch off Info mesaage from chi2 test
@@ -2313,9 +3033,9 @@ int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT)
    
    std::string option = "WW OF UF";
    const char * opt = option.c_str(); 
-   differents |= (h1->Chi2Test(h2, opt) < 1);
-   differents |= (h2->Chi2Test(h1,opt) < 1);         
-   differents |= (bool) equals(h1->Chi2Test(h2,opt), h2->Chi2Test(h1,opt), ERRORLIMIT);
+   differents += (h1->Chi2Test(h2, opt) < 1);
+   differents += (h2->Chi2Test(h1,opt) < 1);         
+   differents += (bool) equals(h1->Chi2Test(h2,opt), h2->Chi2Test(h1,opt), ERRORLIMIT);
    if ( debug )
       cout << "Chi2Test " << h1->Chi2Test(h2, opt) << " " << h2->Chi2Test(h1, opt) 
            << " | " << differents
@@ -2324,7 +3044,7 @@ int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT)
    if (!debug) gErrorIgnoreLevel = precLevel; 
 
    // Mean
-   differents |= (bool) equals(h1->GetMean(1), h2->GetMean(1), ERRORLIMIT);
+   differents += (bool) equals(h1->GetMean(1), h2->GetMean(1), ERRORLIMIT);
    if ( debug )
       cout << "Mean: " << h1->GetMean(1) << " " << h2->GetMean(1) 
            << " | " << fabs( h1->GetMean(1) - h2->GetMean(1) ) 
@@ -2332,7 +3052,7 @@ int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT)
            << endl;
    
    // RMS
-   differents |= (bool) equals( h1->GetRMS(1), h2->GetRMS(1), ERRORLIMIT);
+   differents += (bool) equals( h1->GetRMS(1), h2->GetRMS(1), ERRORLIMIT);
    if ( debug )
       cout << "RMS: " << h1->GetRMS(1) << " " << h2->GetRMS(1) 
            << " | " << fabs( h1->GetRMS(1) - h2->GetRMS(1) ) 
-- 
GitLab