diff --git a/hist/hist/src/TH1.cxx b/hist/hist/src/TH1.cxx
index 387a2bbdd703c321344c2ab342c2aee69da612b1..6fd1985f03f0fa80954ad7f8c48fc37aaa8cbbc3 100644
--- a/hist/hist/src/TH1.cxx
+++ b/hist/hist/src/TH1.cxx
@@ -5292,6 +5292,7 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
    Double_t entries = fEntries;
    Int_t n = TMath::Min(axis->GetNbins(), labels->GetSize());
    std::vector<Int_t> a(n+2);
+   std::vector<Int_t> b(n);
 
    Int_t i,j,k;
    std::vector<Double_t>  cont;
@@ -5312,20 +5313,22 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
       if (GetDimension() == 1) {
          cont.resize(n);
          if (fSumw2.fN) errors.resize(n);
-         for (i=1;i<=n;i++) {
-            cont[i-1] = GetBinContent(i);
-            if (!errors.empty()) errors[i-1] = GetBinError(i);
+         for (i=0; i<n ;i++) {
+            cont[i] = GetBinContent(i+1);
+            if (!errors.empty()) errors[i] = GetBinError(i+1);
+            b[i] = labold->At(i)->GetUniqueID();   // this is the bin corresponding to the label
+            a[i] = i;
          }
          if (sort ==1) TMath::Sort(n,cont.data(),a.data(),kTRUE);  //sort by decreasing values
          else          TMath::Sort(n,cont.data(),a.data(),kFALSE); //sort by increasing values
-         for (i=1;i<=n;i++) {
-            SetBinContent(i,cont[a[i-1]]);
-            if (!errors.empty()) SetBinError(i,errors[a[i-1]]);
+         for (i=0; i<n; i++) {
+            SetBinContent(i+1,cont[b[a[i]]]);
+            if (!errors.empty()) SetBinError(i+1,errors[b[a[i]]]);
          }
-         for (i=1;i<=n;i++) {
-            obj = labold->At(a[i-1]);
+         for (i=0; i<n; i++) {
+            obj = labold->At(a[i]);
             labels->Add(obj);
-            obj->SetUniqueID(i);
+            obj->SetUniqueID(i+1);
          }
       } else if (GetDimension()== 2) {
          std::vector<Double_t> pcont(n+2);
@@ -5371,6 +5374,7 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
       }
    } else {
       //---alphabetic sort
+#if 0
       const UInt_t kUsed = 1<<18;
       TObject *objk=0;
       a[0] = 0;
@@ -5385,6 +5389,9 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
             if (strcmp(label,obj->GetName()) < 0) continue;
             objk = obj;
             a[i] = j;
+            if (gDebug) {
+               Info("LabelsOption","Order labels alhphabetic bin %d goes in %d label %s",j,i);
+            }
             label = obj->GetName();
          }
          if (objk) {
@@ -5398,45 +5405,101 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
          if (!obj) continue;
          obj->ResetBit(kUsed);
       }
+#endif
+      // sort labels using vector of strings and TMath::Sort
+      // I need to array because labels order in list is not necessary that of the bins
+      std::vector<std::string> vecLabels(n);
+      for (i=0; i < n; i++) {
+         vecLabels[i] = labold->At(i)->GetName();
+         b[i] = labold->At(i)->GetUniqueID();   // this is the bin corresponding to the label
+         a[i] = i;
+      }
+      // sort in ascending order for strings
+      TMath::Sort(n, vecLabels.data(), a.data(), kFALSE);
+      // set the new labels
+      for (i=0; i < n; i++) {
+         TObject * labelObj = labold->At(a[i]);
+         labels->Add( labold->At(a[i]) );
+         // set the corresponding bin. NB bin starts from 1
+         labelObj->SetUniqueID(i+1);
+         if (gDebug) std::cout << "bin " << i+1 << " setting new labels for axis " << labold->At(a[i])->GetName() << " from " << b[a[i]] << std::endl;
+      }
 
       if (GetDimension() == 1) {
          cont.resize(n+2);
          if (fSumw2.fN) errors.resize(n+2);
-         for (i=1;i<=n;i++) {
-            cont[i] = GetBinContent(a[i]);
-            if (!errors.empty()) errors[i] = GetBinError(a[i]);
+         for (i=0 ;i < n; i++) {
+            cont[i] = GetBinContent(b[a[i]]);
+            if (!errors.empty()) errors[i] = GetBinError(b[a[i]]);
          }
-         for (i=1;i<=n;i++) {
-            SetBinContent(i,cont[i]);
-            if (!errors.empty()) SetBinError(i,errors[i]);
+         for (i = 0; i < n ; i++) {
+            SetBinContent(i+1,cont[i]);
+            if (!errors.empty()) SetBinError(i+1,errors[i]);
          }
       } else if (GetDimension()== 2) {
          Int_t nx = fXaxis.GetNbins()+2;
          Int_t ny = fYaxis.GetNbins()+2;
          cont.resize(nx*ny);
          if (fSumw2.fN) errors.resize(nx*ny);
-         for (i=0;i<nx;i++) {
-            for (j=0;j<ny;j++) {
-               cont[i+nx*j] = GetBinContent(i,j);
-               if (!errors.empty()) errors[i+nx*j] = GetBinError(i,j);
+         // copy old bin contents and then set to new ordered bins
+         // N.B. bin in histograms starts from 1, but in y we consider under/overflows
+         for (i = 0 ; i < nx; i++) {
+            for (j = 0; j < ny; j++) {  // ny is nbins+2
+               cont[i + nx*j] = GetBinContent(i,j);
+               if (!errors.empty()) errors[i + nx*j] = GetBinError(i,j);
             }
          }
          if (axis == GetXaxis()) {
-            for (i=1;i<=n;i++) {
-               for (j=0;j<ny;j++) {
-                  SetBinContent(i,j,cont[a[i]+nx*j]);
-                  if (!errors.empty()) SetBinError(i,j,errors[a[i]+nx*j]);
+            for (i = 0; i < n; i++) {
+               for (j = 0; j < ny; j++) {
+                  SetBinContent(i + 1, j, cont[b[a[i]] + nx * j]);
+                  if (!errors.empty())
+                     SetBinError(i + 1, j, errors[b[a[i]] + nx * j]);
                }
             }
          } else {
-            for (i=0;i<nx;i++) {
-               for (j=1;j<=n;j++) {
-                  SetBinContent(i,j,cont[i+nx*a[j]]);
-                  if (!errors.empty()) SetBinError(i,j,errors[i+nx*a[j]]);
+            for (i = 0; i < nx; i++) {
+               for (j = 0; j < n; j++) {
+                  SetBinContent(i, j + 1, cont[i + nx * b[a[j]]]);
+                  if (!errors.empty())
+                     SetBinError(i, j + 1, errors[i + nx * b[a[j]]]);
                }
             }
          }
+#if 0
+         if (axis == GetXaxis()) {
+            // copy old bin contents and then set to new ordered bins
+            // N.B. bin in histograms starts from 1, but in y we consider under/overflows
+            for (i = 0 ; i < n; i++) {
+               for (j = 0; j < ny; j++) {  // ny is nbins+2
+                  cont[i + n*j] = GetBinContent(b[a[i]],j);
+                  if (!errors.empty()) errors[i + n*j] = GetBinError( b[a[i]],j);
+               }
+            }
+            for (i = 0; i < n; i++) {
+               for (j = 0; j< ny; j++) {
+                  SetBinContent(i+1, j, cont[i + n*j]);
+                  if (!errors.empty()) SetBinError(i+1 , j, errors[i + n*j]);
+               }
+            }
+         } else {
+             for (i = 0 ; i < nx; i++) {
+               for (j = 0; j < n; j++) {
+                  cont[j + n*i] = GetBinContent(i,b[a[j]]);
+                  if (!errors.empty()) errors[j + n*i] = GetBinError( i, b[a[j]]);
+               }
+            }
+            for (i = 0; i < nx; i++) {
+               for (j = 0; j< n; j++) {
+                  SetBinContent(i, j+1, cont[j + n*i]);
+                  if (!errors.empty()) SetBinError(i, j+1, errors[j + n*i]);
+               }
+            }
+         }
+#endif
+
       } else {
+         // case of 3D (needs to be fixed!)
          Int_t nx = fXaxis.GetNbins()+2;
          Int_t ny = fYaxis.GetNbins()+2;
          Int_t nz = fZaxis.GetNbins()+2;
@@ -5452,11 +5515,11 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
          }
          if (axis == GetXaxis()) {
             // labels on x axis
-            for (i=1;i<=n;i++) {
+            for (i=0; i<n; i++) {  // for x we lool only on bins with the labels
                for (j=0;j<ny;j++) {
                   for (k=0;k<nz;k++) {
-                     SetBinContent(i,j,k,cont[a[i]+nx*(j+ny*k)]);
-                     if (!errors.empty()) SetBinError(i,j,k,errors[a[i]+nx*(j+ny*k)]);
+                     SetBinContent(i+1,j,k,cont[b[a[i]]+nx*(j+ny*k)]);
+                     if (!errors.empty()) SetBinError(i+1,j,k,errors[b[a[i]]+nx*(j+ny*k)]);
                   }
                }
             }
@@ -5464,10 +5527,10 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
          else if (axis == GetYaxis()) {
             // labels on y axis
             for (i=0;i<nx;i++) {
-               for (j=1;j<=n;j++) {
+               for (j=0; j<n; j++) {
                   for (k=0;k<nz;k++) {
-                     SetBinContent(i,j,k,cont[i+nx*(a[j]+ny*k)]);
-                     if (!errors.empty()) SetBinError(i,j,k,errors[i+nx*(a[j]+ny*k)]);
+                     SetBinContent(i,j+1,k,cont[i+nx*(b[a[j]]+ny*k)]);
+                     if (!errors.empty()) SetBinError(i,j+1,k,errors[i+nx*(b[a[j]]+ny*k)]);
                   }
                }
             }
@@ -5476,15 +5539,17 @@ void TH1::LabelsOption(Option_t *option, Option_t *ax)
             // labels on z axis
             for (i=0;i<nx;i++) {
                for (j=0;j<ny;j++) {
-                  for (k=1;k<=n;k++) {
-                     SetBinContent(i,j,k,cont[i+nx*(j+ny*a[k])]);
-                     if (!errors.empty()) SetBinError(i,j,k,errors[i+nx*(j+ny*a[k])]);
+                  for (k=0; k<n; k++) {
+                     SetBinContent(i,j,k+1,cont[i+nx*(j+ny*b[a[k]])]);
+                     if (!errors.empty()) SetBinError(i,j,k+1,errors[i+nx*(j+ny*b[a[k]])]);
                   }
                }
             }
          }
       }
    }
+   // need to reset stattistics after sorting
+   ResetStats();
    fEntries = entries;
    delete labold;
 }
@@ -7373,7 +7438,7 @@ void TH1::GetStats(Double_t *stats) const
    Double_t x;
    // identify the case of labels with extension of axis range
    // in this case the statistics in x does not make any sense
-   Bool_t labelHist =  ((const_cast<TAxis&>(fXaxis)).GetLabels() && CanExtendAllAxes() );
+   Bool_t labelHist =  ((const_cast<TAxis&>(fXaxis)).GetLabels() && fXaxis.CanExtend() );
    // fTsumw == 0 && fEntries > 0 is a special case when uses SetBinContent or calls ResetStats before
    if ((fTsumw == 0 && fEntries > 0) || ( fXaxis.TestBit(TAxis::kAxisRange) && !labelHist) ) {
       for (bin=0;bin<4;bin++) stats[bin] = 0;
diff --git a/hist/hist/src/TH2.cxx b/hist/hist/src/TH2.cxx
index ada1155483975f5c98e2f4576e3c3ad5c50c1e4f..4fcf50408c59be3caf8c7ac8df9d025594e9235f 100644
--- a/hist/hist/src/TH2.cxx
+++ b/hist/hist/src/TH2.cxx
@@ -717,7 +717,7 @@ void TH2::DoFitSlices(bool onX,
    if (lastbin < 0 || lastbin > nbins + 1) lastbin = nbins + 1;
    if (lastbin < firstbin) {firstbin = 0; lastbin = nbins + 1;}
 
-   
+
    TString opt = option;
    TString proj_opt = "e";
    Int_t i1 = opt.Index("[");
@@ -761,11 +761,11 @@ void TH2::DoFitSlices(bool onX,
    const TArrayD *bins = outerAxis.GetXbins();
    // outer axis boudaries used for creating reported histograms are different
    // than the limits used in the projection loop (firstbin,lastbin)
-   Int_t firstOutBin = outerAxis.TestBit(TAxis::kAxisRange) ? std::max(firstbin,1) : 1; 
+   Int_t firstOutBin = outerAxis.TestBit(TAxis::kAxisRange) ? std::max(firstbin,1) : 1;
    Int_t lastOutBin = outerAxis.TestBit(TAxis::kAxisRange) ?  std::min(lastbin,outerAxis.GetNbins() ) : outerAxis.GetNbins();
    Int_t nOutBins = lastOutBin-firstOutBin+1;
    // merge bins if use nstep > 1 and fixed bins
-   if (bins->fN == 0) nOutBins /= nstep;  
+   if (bins->fN == 0) nOutBins /= nstep;
    for (ipar=0;ipar<npar;ipar++) {
       snprintf(name,2000,"%s_%d",GetName(),ipar);
       snprintf(title,2000,"Fitted value of par[%d]=%s",ipar,f1->GetParName(ipar));
@@ -831,7 +831,7 @@ void TH2::DoFitSlices(bool onX,
       }
       // don't need to delete hp. If histogram has the same name it is re-used in TH2::Projection
    }
-   delete hp; 
+   delete hp;
    delete [] parsave;
    delete [] name;
    delete [] title;
@@ -1148,10 +1148,14 @@ void TH2::GetStats(Double_t *stats) const
             if (lastBinY ==  fYaxis.GetNbins() ) lastBinY += 1;
          }
       }
+      // check for labels axis . In that case corresponsing statistics do not make sense and it is set to zero
+      Bool_t labelXaxis =  ((const_cast<TAxis&>(fXaxis)).GetLabels() && fXaxis.CanExtend() );
+      Bool_t labelYaxis =  ((const_cast<TAxis&>(fYaxis)).GetLabels() && fYaxis.CanExtend() );
+
       for (Int_t biny = firstBinY; biny <= lastBinY; ++biny) {
-         Double_t y = fYaxis.GetBinCenter(biny);
+         Double_t y = (!labelYaxis) ? fYaxis.GetBinCenter(biny) : 0;
          for (Int_t binx = firstBinX; binx <= lastBinX; ++binx) {
-            Double_t x = fXaxis.GetBinCenter(binx);
+            Double_t x = (!labelXaxis) ? fXaxis.GetBinCenter(binx) : 0;
             //w   = TMath::Abs(GetBinContent(bin));
             Int_t bin = GetBin(binx,biny);
             Double_t w = RetrieveBinContent(bin);
diff --git a/hist/hist/src/TH3.cxx b/hist/hist/src/TH3.cxx
index c0a7e94ebae6931d9ec0014f5348d3a83c00832c..ebfaac09520bbc598f8e910717ab933c24883384 100644
--- a/hist/hist/src/TH3.cxx
+++ b/hist/hist/src/TH3.cxx
@@ -1163,13 +1163,19 @@ void TH3::GetStats(Double_t *stats) const
             if (lastBinZ ==  fZaxis.GetNbins() ) lastBinZ += 1;
          }
       }
+
+      // check for labels axis . In that case corresponsing statistics do not make sense and it is set to zero
+      Bool_t labelXaxis =  ((const_cast<TAxis&>(fXaxis)).GetLabels() && fXaxis.CanExtend() );
+      Bool_t labelYaxis =  ((const_cast<TAxis&>(fYaxis)).GetLabels() && fYaxis.CanExtend() );
+      Bool_t labelZaxis =  ((const_cast<TAxis&>(fZaxis)).GetLabels() && fZaxis.CanExtend() );
+
       for (binz = firstBinZ; binz <= lastBinZ; binz++) {
-         z = fZaxis.GetBinCenter(binz);
+         z = (!labelZaxis) ? fZaxis.GetBinCenter(binz) : 0;
          for (biny = firstBinY; biny <= lastBinY; biny++) {
-            y = fYaxis.GetBinCenter(biny);
+            y = (!labelYaxis) ? fYaxis.GetBinCenter(biny) : 0;
             for (binx = firstBinX; binx <= lastBinX; binx++) {
                bin = GetBin(binx,biny,binz);
-               x   = fXaxis.GetBinCenter(binx);
+               x = (!labelXaxis) ? fXaxis.GetBinCenter(binx) : 0;
                //w   = TMath::Abs(GetBinContent(bin));
                w   = RetrieveBinContent(bin);
                err = TMath::Abs(GetBinError(bin));
diff --git a/test/stressHistogram.cxx b/test/stressHistogram.cxx
index 8212147085cfb040f608974665f12ae58f1c2cbe..0e66af2905703d5073f520ab693c3e46605eb9a8 100644
--- a/test/stressHistogram.cxx
+++ b/test/stressHistogram.cxx
@@ -94,6 +94,7 @@
 
 #include "TROOT.h"
 #include <algorithm>
+#include <random>
 #include <cassert>
 
 using namespace std;
@@ -158,6 +159,7 @@ int equals(const char* msg, TH3D* h1, TH3D* h2, int options = 0, double ERRORLIM
 int equals(const char* msg, THnBase* h1, THnBase* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
 int equals(const char* msg, THnBase* h1, TH1* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
 int equals(Double_t n1, Double_t n2, double ERRORLIMIT = defaultErrorLimit);
+int equals(const char * s1, const char * s2);  // for comparing names (e.g. axis labels)
 int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT = defaultErrorLimit);
 std::ostream& operator<<(std::ostream& out, TH1D* h);
 // old stresHistOpts.cxx file
@@ -169,16 +171,16 @@ bool testAdd1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH1D* h1 = new TH1D("t1D1-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("t1D1-h2", "h2-Title", numberOfBins, minRange, maxRange);
-   TH1D* h3 = new TH1D("t1D1-h3", "h3=c1*h1+c2*h2", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("t1D1_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("t1D1_h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("t1D1_h3", "h3=c1*h1+c2*h2", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();h2->Sumw2();h3->Sumw2();
 
    FillHistograms(h1, h3, 1.0, c1);
    FillHistograms(h2, h3, 1.0, c2);
 
-   TH1D* h4 = new TH1D("t1D1-h4", "h4=c1*h1+h2*c2", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("t1D1_h4", "h4=c1*h1+h2*c2", numberOfBins, minRange, maxRange);
    h4->Add(h1, h2, c1, c2);
 
    bool ret = equals("Add1D1", h3, h4, cmpOptStats, 1E-13);
@@ -284,8 +286,8 @@ bool testAdd3()
 
    Double_t c1 = r.Rndm();
 
-   TH1D* h1 = new TH1D("t1D1-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("t1D1-h2", "h2=c1*h1+c2*h2", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("t1D1_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("t1D1_h2", "h2=c1*h1+c2*h2", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();h2->Sumw2();
 
@@ -296,7 +298,7 @@ bool testAdd3()
    }
 
 
-   TH1D* h3 = new TH1D("t1D1-h3", "h3=c1*h1", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("t1D1_h3", "h3=c1*h1", numberOfBins, minRange, maxRange);
    h3->Add(h1, h1, c1, -1);
 
    // TH1::Add will reset the stats in this case so we need to do for the reference histogram
@@ -327,7 +329,7 @@ bool testAddVar1()
    FillHistograms(h1, h3, 1.0, c1);
    FillHistograms(h2, h3, 1.0, c2);
 
-   TH1D* h4 = new TH1D("t1D1-h4", "h4=c1*h1+h2*c2", numberOfBins, v);
+   TH1D* h4 = new TH1D("t1D1_h4", "h4=c1*h1+h2*c2", numberOfBins, v);
    h4->Add(h1, h2, c1, c2);
 
    bool ret = equals("AddVar1D1", h3, h4, cmpOptStats, 1E-13);
@@ -427,8 +429,8 @@ bool testAddVar3()
 
    Double_t c1 = r.Rndm();
 
-   TH1D* h1 = new TH1D("t1D1-h1", "h1-Title", numberOfBins, v);
-   TH1D* h2 = new TH1D("t1D1-h2", "h2=c1*h1+c2*h2", numberOfBins, v);
+   TH1D* h1 = new TH1D("t1D1_h1", "h1-Title", numberOfBins, v);
+   TH1D* h2 = new TH1D("t1D1_h2", "h2=c1*h1+c2*h2", numberOfBins, v);
 
    h1->Sumw2();h2->Sumw2();
 
@@ -438,7 +440,7 @@ bool testAddVar3()
       h2->Fill(value, c1 / h1->GetBinWidth( h1->FindBin(value) ) );
    }
 
-   TH1D* h3 = new TH1D("t1D1-h3", "h3=c1*h1", numberOfBins, v);
+   TH1D* h3 = new TH1D("t1D1_h3", "h3=c1*h1", numberOfBins, v);
    h3->Add(h1, h1, c1, -1);
 
    // TH1::Add will reset the stats in this case so we need to do for the reference histogram
@@ -457,10 +459,10 @@ bool testAdd2D3()
 
    Double_t c1 = r.Rndm();
 
-   TH2D* h1 = new TH2D("t1D1-h1", "h1-Title",
+   TH2D* h1 = new TH2D("t1D1_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins+2, minRange, maxRange);
-   TH2D* h2 = new TH2D("t1D1-h2", "h2=c1*h1+c2*h2",
+   TH2D* h2 = new TH2D("t1D1_h2", "h2=c1*h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins+2, minRange, maxRange);
 
@@ -476,7 +478,7 @@ bool testAdd2D3()
       h2->Fill(x, y, c1 / area);
    }
 
-   TH2D* h3 = new TH2D("t1D1-h3", "h3=c1*h1",
+   TH2D* h3 = new TH2D("t1D1_h3", "h3=c1*h1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins+2, minRange, maxRange);
    h3->Add(h1, h1, c1, -1);
@@ -496,11 +498,11 @@ bool testAdd3D3()
 
    Double_t c1 = r.Rndm();
 
-   TH3D* h1 = new TH3D("t1D1-h1", "h1-Title",
+   TH3D* h1 = new TH3D("t1D1_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins+1, minRange, maxRange,
                        numberOfBins+2, minRange, maxRange);
-   TH3D* h2 = new TH3D("t1D1-h2", "h2=c1*h1+c2*h2",
+   TH3D* h2 = new TH3D("t1D1_h2", "h2=c1*h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins+1, minRange, maxRange,
                        numberOfBins+2, minRange, maxRange);
@@ -521,7 +523,7 @@ bool testAdd3D3()
       h2->Fill(x, y, z, c1 / area);
    }
 
-   TH3D* h3 = new TH3D("t1D1-h3", "h3=c1*h1",
+   TH3D* h3 = new TH3D("t1D1_h3", "h3=c1*h1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins+1, minRange, maxRange,
                        numberOfBins+2, minRange, maxRange);
@@ -543,15 +545,15 @@ bool testAdd2D1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH2D* h1 = new TH2D("t2D1-h1", "h1",
+   TH2D* h1 = new TH2D("t2D1_h1", "h1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH2D* h2 = new TH2D("t2D1-h2", "h2",
+   TH2D* h2 = new TH2D("t2D1_h2", "h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH2D* h3 = new TH2D("t2D1-h3", "h3=c1*h1+c2*h2",
+   TH2D* h3 = new TH2D("t2D1_h3", "h3=c1*h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -571,7 +573,7 @@ bool testAdd2D1()
       h3->Fill(x, y, c2);
    }
 
-   TH2D* h4 = new TH2D("t2D1-h4", "h4=c1*h1+c2*h2",
+   TH2D* h4 = new TH2D("t2D1_h4", "h4=c1*h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    h4->Add(h1, h2, c1, c2);
@@ -634,15 +636,15 @@ bool testAdd2D2()
 
    Double_t c2 = r.Rndm();
 
-   TH2D* h1 = new TH2D("t2D2-h1", "h1",
+   TH2D* h1 = new TH2D("t2D2_h1", "h1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH2D* h2 = new TH2D("t2D2-h2", "h2",
+   TH2D* h2 = new TH2D("t2D2_h2", "h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH2D* h3 = new TH2D("t2D2-h3", "h3=h1+c2*h2",
+   TH2D* h3 = new TH2D("t2D2_h3", "h3=h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -717,17 +719,17 @@ bool testAdd3D1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH3D* h1 = new TH3D("t3D1-h1", "h1",
+   TH3D* h1 = new TH3D("t3D1_h1", "h1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH3D* h2 = new TH3D("t3D1-h2", "h2",
+   TH3D* h2 = new TH3D("t3D1_h2", "h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH3D* h3 = new TH3D("t3D1-h3", "h3=c1*h1+c2*h2",
+   TH3D* h3 = new TH3D("t3D1_h3", "h3=c1*h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -750,7 +752,7 @@ bool testAdd3D1()
       h3->Fill(x, y, z, c2);
    }
 
-   TH3D* h4 = new TH3D("t3D1-h4", "h4=c1*h1+c2*h2",
+   TH3D* h4 = new TH3D("t3D1_h4", "h4=c1*h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -820,17 +822,17 @@ bool testAdd3D2()
 
    Double_t c2 = r.Rndm();
 
-   TH3D* h1 = new TH3D("t3D2-h1", "h1",
+   TH3D* h1 = new TH3D("t3D2_h1", "h1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH3D* h2 = new TH3D("t3D2-h2", "h2",
+   TH3D* h2 = new TH3D("t3D2_h2", "h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
-   TH3D* h3 = new TH3D("t3D2-h3", "h3=h1+c2*h2",
+   TH3D* h3 = new TH3D("t3D2_h3", "h3=h1+c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -957,9 +959,9 @@ bool testMul1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH1D* h1 = new TH1D("m1D1-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("m1D1-h2", "h2-Title", numberOfBins, minRange, maxRange);
-   TH1D* h3 = new TH1D("m1D1-h3", "h3=c1*h1*c2*h2", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("m1D1_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("m1D1_h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("m1D1_h3", "h3=c1*h1*c2*h2", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();h2->Sumw2();h3->Sumw2();
 
@@ -989,7 +991,7 @@ bool testMul1()
       h3->SetBinContent(bin, h3->GetBinContent(bin) / 2 );
    }
 
-   TH1D* h4 = new TH1D("m1D1-h4", "h4=h1*h2", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("m1D1_h4", "h4=h1*h2", numberOfBins, minRange, maxRange);
    h4->Multiply(h1, h2, c1, c2);
 
    bool ret = equals("Multiply1D1", h3, h4, cmpOptStats  , 1E-14);
@@ -1009,9 +1011,9 @@ bool testMulVar1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH1D* h1 = new TH1D("m1D1-h1", "h1-Title", numberOfBins, v);
-   TH1D* h2 = new TH1D("m1D1-h2", "h2-Title", numberOfBins, v);
-   TH1D* h3 = new TH1D("m1D1-h3", "h3=c1*h1*c2*h2", numberOfBins, v);
+   TH1D* h1 = new TH1D("m1D1_h1", "h1-Title", numberOfBins, v);
+   TH1D* h2 = new TH1D("m1D1_h2", "h2-Title", numberOfBins, v);
+   TH1D* h3 = new TH1D("m1D1_h3", "h3=c1*h1*c2*h2", numberOfBins, v);
 
    h1->Sumw2();h2->Sumw2();h3->Sumw2();
 
@@ -1041,7 +1043,7 @@ bool testMulVar1()
       h3->SetBinContent(bin, h3->GetBinContent(bin) / 2 );
    }
 
-   TH1D* h4 = new TH1D("m1D1-h4", "h4=h1*h2", numberOfBins, v);
+   TH1D* h4 = new TH1D("m1D1_h4", "h4=h1*h2", numberOfBins, v);
    h4->Multiply(h1, h2, c1, c2);
 
    bool ret = equals("MultiVar1D1", h3, h4, cmpOptStats, 1E-14);
@@ -1055,9 +1057,9 @@ bool testMul2()
 {
    // Tests the second Multiply method for 1D Histograms
 
-   TH1D* h1 = new TH1D("m1D2-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("m1D2-h2", "h2-Title", numberOfBins, minRange, maxRange);
-   TH1D* h3 = new TH1D("m1D2-h3", "h3=h1*h2", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("m1D2_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("m1D2_h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("m1D2_h3", "h3=h1*h2", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();h2->Sumw2();h3->Sumw2();
 
@@ -1100,9 +1102,9 @@ bool testMulVar2()
    Double_t v[numberOfBins+1];
    FillVariableRange(v);
 
-   TH1D* h1 = new TH1D("m1D2-h1", "h1-Title", numberOfBins, v);
-   TH1D* h2 = new TH1D("m1D2-h2", "h2-Title", numberOfBins, v);
-   TH1D* h3 = new TH1D("m1D2-h3", "h3=h1*h2", numberOfBins, v);
+   TH1D* h1 = new TH1D("m1D2_h1", "h1-Title", numberOfBins, v);
+   TH1D* h2 = new TH1D("m1D2_h2", "h2-Title", numberOfBins, v);
+   TH1D* h3 = new TH1D("m1D2_h3", "h3=h1*h2", numberOfBins, v);
 
    h1->Sumw2();h2->Sumw2();h3->Sumw2();
 
@@ -1145,13 +1147,13 @@ bool testMul2D1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH2D* h1 = new TH2D("m2D1-h1", "h1-Title",
+   TH2D* h1 = new TH2D("m2D1_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("m2D1-h2", "h2-Title",
+   TH2D* h2 = new TH2D("m2D1_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h3 = new TH2D("m2D1-h3", "h3=c1*h1*c2*h2",
+   TH2D* h3 = new TH2D("m2D1_h3", "h3=c1*h1*c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -1190,7 +1192,7 @@ bool testMul2D1()
       }
    }
 
-   TH2D* h4 = new TH2D("m2D1-h4", "h4=h1*h2",
+   TH2D* h4 = new TH2D("m2D1_h4", "h4=h1*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    h4->Multiply(h1, h2, c1, c2);
@@ -1206,13 +1208,13 @@ bool testMul2D2()
 {
    // Tests the second Multiply method for 2D Histograms
 
-   TH2D* h1 = new TH2D("m2D2-h1", "h1-Title",
+   TH2D* h1 = new TH2D("m2D2_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("m2D2-h2", "h2-Title",
+   TH2D* h2 = new TH2D("m2D2_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h3 = new TH2D("m2D2-h3", "h3=h1*h2",
+   TH2D* h3 = new TH2D("m2D2_h3", "h3=h1*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -1266,15 +1268,15 @@ bool testMul3D1()
    Double_t c1 = r.Rndm();
    Double_t c2 = r.Rndm();
 
-   TH3D* h1 = new TH3D("m3D1-h1", "h1-Title",
+   TH3D* h1 = new TH3D("m3D1_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("m3D1-h2", "h2-Title",
+   TH3D* h2 = new TH3D("m3D1_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("m3D1-h3", "h3=c1*h1*c2*h2",
+   TH3D* h3 = new TH3D("m3D1_h3", "h3=c1*h1*c2*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -1321,7 +1323,7 @@ bool testMul3D1()
       }
    }
 
-   TH3D* h4 = new TH3D("m3D1-h4", "h4=h1*h2",
+   TH3D* h4 = new TH3D("m3D1_h4", "h4=h1*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -1338,15 +1340,15 @@ bool testMul3D2()
 {
    // Tests the second Multiply method for 3D Histograms
 
-   TH3D* h1 = new TH3D("m3D2-h1", "h1-Title",
+   TH3D* h1 = new TH3D("m3D2_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("m3D2-h2", "h2-Title",
+   TH3D* h2 = new TH3D("m3D2_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("m3D2-h3", "h3=h1*h2",
+   TH3D* h3 = new TH3D("m3D2_h3", "h3=h1*h2",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -1475,8 +1477,8 @@ bool testMulF1D()
 {
    Double_t c1 = r.Rndm();
 
-   TH1D* h1 = new TH1D("mf1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("mf1D-h2", "h2=h1*c1*f1", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("mf1D_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("mf1D_h2", "h2=h1*c1*f1", numberOfBins, minRange, maxRange);
 
    TF1* f = new TF1("sin", "sin(x)", minRange - 2, maxRange + 2);
 
@@ -1504,8 +1506,8 @@ bool testMulF1D2()
 {
    Double_t c1 = r.Rndm();
 
-   TH1D* h1 = new TH1D("mf1D2-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("mf1D2-h2", "h2=h1*c1*f1", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("mf1D2_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("mf1D2_h2", "h2=h1*c1*f1", numberOfBins, minRange, maxRange);
 
    TF2* f = new TF2("sin2", "sin(x)*cos(y)",
                     minRange - 2, maxRange + 2,
@@ -1536,10 +1538,10 @@ bool testMulF2D()
 {
    Double_t c1 = r.Rndm();
 
-   TH2D* h1 = new TH2D("mf2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("mf2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
-   TH2D* h2 = new TH2D("mf2D-h2", "h2=h1*c1*f1",
+   TH2D* h2 = new TH2D("mf2D_h2", "h2=h1*c1*f1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
 
@@ -1570,10 +1572,10 @@ bool testMulF2D2()
 {
    Double_t c1 = r.Rndm();
 
-   TH2D* h1 = new TH2D("mf2D2-h1", "h1-Title",
+   TH2D* h1 = new TH2D("mf2D2_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
-   TH2D* h2 = new TH2D("mf2D2-h2", "h2=h1*c1*f1",
+   TH2D* h2 = new TH2D("mf2D2_h2", "h2=h1*c1*f1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
 
@@ -1608,11 +1610,11 @@ bool testMulF3D()
 {
    Double_t c1 = r.Rndm();
 
-   TH3D* h1 = new TH3D("mf3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("mf3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
-   TH3D* h2 = new TH3D("mf3D-h2", "h2=h1*c1*f1",
+   TH3D* h2 = new TH3D("mf3D_h2", "h2=h1*c1*f1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
@@ -1645,11 +1647,11 @@ bool testMulF3D2()
 {
    Double_t c1 = r.Rndm();
 
-   TH3D* h1 = new TH3D("mf3D2-h1", "h1-Title",
+   TH3D* h1 = new TH3D("mf3D2_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
-   TH3D* h2 = new TH3D("mf3D2-h2", "h2=h1*c1*f1",
+   TH3D* h2 = new TH3D("mf3D2_h2", "h2=h1*c1*f1",
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange,
                        numberOfBins, minRange, maxRange);
@@ -1769,8 +1771,8 @@ bool testDivide1()
    Double_t c1 = r.Rndm() + 1;
    Double_t c2 = r.Rndm() + 1;
 
-   TH1D* h1 = new TH1D("d1D1-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("d1D1-h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("d1D1_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("d1D1_h2", "h2-Title", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();h2->Sumw2();
 
@@ -1789,10 +1791,10 @@ bool testDivide1()
       if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);
 
 
-   TH1D* h3 = new TH1D("d1D1-h3", "h3=(c1*h1)/(c2*h2)", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("d1D1_h3", "h3=(c1*h1)/(c2*h2)", numberOfBins, minRange, maxRange);
    h3->Divide(h1, h2, c1, c2);
 
-   TH1D* h4 = new TH1D("d1D1-h4", "h4=h3*h2)", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("d1D1_h4", "h4=h3*h2)", numberOfBins, minRange, maxRange);
    h4->Multiply(h2, h3, c2/c1, 1);
    for ( Int_t bin = 0; bin <= h4->GetNbinsX() + 1; ++bin ) {
       Double_t error = h4->GetBinError(bin) * h4->GetBinError(bin);
@@ -1819,8 +1821,8 @@ bool testDivideVar1()
    Double_t c1 = r.Rndm() + 1;
    Double_t c2 = r.Rndm() + 1;
 
-   TH1D* h1 = new TH1D("d1D1-h1", "h1-Title", numberOfBins, v);
-   TH1D* h2 = new TH1D("d1D1-h2", "h2-Title", numberOfBins, v);
+   TH1D* h1 = new TH1D("d1D1_h1", "h1-Title", numberOfBins, v);
+   TH1D* h2 = new TH1D("d1D1_h2", "h2-Title", numberOfBins, v);
 
    h1->Sumw2();h2->Sumw2();
 
@@ -1839,10 +1841,10 @@ bool testDivideVar1()
       if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);
 
 
-   TH1D* h3 = new TH1D("d1D1-h3", "h3=(c1*h1)/(c2*h2)", numberOfBins, v);
+   TH1D* h3 = new TH1D("d1D1_h3", "h3=(c1*h1)/(c2*h2)", numberOfBins, v);
    h3->Divide(h1, h2, c1, c2);
 
-   TH1D* h4 = new TH1D("d1D1-h4", "h4=h3*h2)", numberOfBins, v);
+   TH1D* h4 = new TH1D("d1D1_h4", "h4=h3*h2)", numberOfBins, v);
    h4->Multiply(h2, h3, c2/c1, 1);
    for ( Int_t bin = 0; bin <= h4->GetNbinsX() + 1; ++bin ) {
       Double_t error = h4->GetBinError(bin) * h4->GetBinError(bin);
@@ -1901,8 +1903,8 @@ bool testDivide2()
 {
    // Tests the second Divide method for 1D Histograms
 
-   TH1D* h1 = new TH1D("d1D2-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("d1D2-h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("d1D2_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("d1D2_h2", "h2-Title", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();h2->Sumw2();
 
@@ -1923,7 +1925,7 @@ bool testDivide2()
    TH1D* h3 = static_cast<TH1D*>( h1->Clone() );
    h3->Divide(h2);
 
-   TH1D* h4 = new TH1D("d1D2-h4", "h4=h3*h2)", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("d1D2_h4", "h4=h3*h2)", numberOfBins, minRange, maxRange);
    h4->Multiply(h2, h3, 1.0, 1.0);
    for ( Int_t bin = 0; bin <= h4->GetNbinsX() + 1; ++bin ) {
       Double_t error = h4->GetBinError(bin) * h4->GetBinError(bin);
@@ -1948,8 +1950,8 @@ bool testDivideVar2()
    Double_t v[numberOfBins+1];
    FillVariableRange(v);
 
-   TH1D* h1 = new TH1D("d1D2-h1", "h1-Title", numberOfBins, v);
-   TH1D* h2 = new TH1D("d1D2-h2", "h2-Title", numberOfBins, v);
+   TH1D* h1 = new TH1D("d1D2_h1", "h1-Title", numberOfBins, v);
+   TH1D* h2 = new TH1D("d1D2_h2", "h2-Title", numberOfBins, v);
 
    h1->Sumw2();h2->Sumw2();
 
@@ -1970,7 +1972,7 @@ bool testDivideVar2()
    TH1D* h3 = static_cast<TH1D*>( h1->Clone() );
    h3->Divide(h2);
 
-   TH1D* h4 = new TH1D("d1D2-h4", "h4=h3*h2)", numberOfBins, v);
+   TH1D* h4 = new TH1D("d1D2_h4", "h4=h3*h2)", numberOfBins, v);
    h4->Multiply(h2, h3, 1.0, 1.0);
    for ( Int_t bin = 0; bin <= h4->GetNbinsX() + 1; ++bin ) {
       Double_t error = h4->GetBinError(bin) * h4->GetBinError(bin);
@@ -1995,10 +1997,10 @@ bool testDivide2D1()
    Double_t c1 = r.Rndm() + 1;
    Double_t c2 = r.Rndm() + 1;
 
-   TH2D* h1 = new TH2D("d2D1-h1", "h1-Title",
+   TH2D* h1 = new TH2D("d2D1_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("d2D1-h2", "h2-Title",
+   TH2D* h2 = new TH2D("d2D1_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -2020,12 +2022,12 @@ bool testDivide2D1()
    for (int i = 0; i < h2->GetSize(); ++i)
       if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);
 
-   TH2D* h3 = new TH2D("d2D1-h3", "h3=(c1*h1)/(c2*h2)",
+   TH2D* h3 = new TH2D("d2D1_h3", "h3=(c1*h1)/(c2*h2)",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    h3->Divide(h1, h2, c1, c2);
 
-   TH2D* h4 = new TH2D("d2D1-h4", "h4=h3*h2)",
+   TH2D* h4 = new TH2D("d2D1_h4", "h4=h3*h2)",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    h4->Multiply(h2, h3, c2/c1, 1);
@@ -2051,10 +2053,10 @@ bool testDivide2D2()
 {
    // Tests the second Divide method for 2D Histograms
 
-   TH2D* h1 = new TH2D("d2D2-h1", "h1-Title",
+   TH2D* h1 = new TH2D("d2D2_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("d2D2-h2", "h2-Title",
+   TH2D* h2 = new TH2D("d2D2_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -2079,7 +2081,7 @@ bool testDivide2D2()
    TH2D* h3 = static_cast<TH2D*>( h1->Clone() );
    h3->Divide(h2);
 
-   TH2D* h4 = new TH2D("d2D2-h4", "h4=h3*h2)",
+   TH2D* h4 = new TH2D("d2D2_h4", "h4=h3*h2)",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    h4->Multiply(h2, h3, 1.0, 1.0);
@@ -2108,11 +2110,11 @@ bool testDivide3D1()
    Double_t c1 = r.Rndm() + 1;
    Double_t c2 = r.Rndm() + 1;
 
-   TH3D* h1 = new TH3D("d3D1-h1", "h1-Title",
+   TH3D* h1 = new TH3D("d3D1_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("d3D1-h2", "h2-Title",
+   TH3D* h2 = new TH3D("d3D1_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2137,13 +2139,13 @@ bool testDivide3D1()
    for (int i = 0; i < h2->GetSize(); ++i)
       if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);
 
-   TH3D* h3 = new TH3D("d3D1-h3", "h3=(c1*h1)/(c2*h2)",
+   TH3D* h3 = new TH3D("d3D1_h3", "h3=(c1*h1)/(c2*h2)",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    h3->Divide(h1, h2, c1, c2);
 
-   TH3D* h4 = new TH3D("d3D1-h4", "h4=h3*h2)",
+   TH3D* h4 = new TH3D("d3D1_h4", "h4=h3*h2)",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2174,11 +2176,11 @@ bool testDivide3D2()
 {
    // Tests the second Divide method for 3D Histograms
 
-   TH3D* h1 = new TH3D("d3D2-h1", "h1-Title",
+   TH3D* h1 = new TH3D("d3D2_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("d3D2-h2", "h2-Title",
+   TH3D* h2 = new TH3D("d3D2_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2206,7 +2208,7 @@ bool testDivide3D2()
    TH3D* h3 = static_cast<TH3D*>( h1->Clone() );
    h3->Divide(h2);
 
-   TH3D* h4 = new TH3D("d3D2-h4", "h4=h3*h2)",
+   TH3D* h4 = new TH3D("d3D2_h4", "h4=h3*h2)",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2363,7 +2365,7 @@ bool testAssign1D()
 {
    // Tests the operator=() method for 1D Histograms
 
-   TH1D* h1 = new TH1D("=1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("=1D_h1", "h1-Title", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();
 
@@ -2372,7 +2374,7 @@ bool testAssign1D()
       h1->Fill(value, 1.0);
    }
 
-   TH1D* h2 = new TH1D("=1D-h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("=1D_h2", "h2-Title", numberOfBins, minRange, maxRange);
    *h2 = *h1;
 
    bool ret = equals("Assign Oper Hist '='  1D", h1, h2, cmpOptStats);
@@ -2387,7 +2389,7 @@ bool testAssignVar1D()
    Double_t v[numberOfBins+1];
    FillVariableRange(v);
 
-   TH1D* h1 = new TH1D("=1D-h1", "h1-Title", numberOfBins, v);
+   TH1D* h1 = new TH1D("=1D_h1", "h1-Title", numberOfBins, v);
 
    h1->Sumw2();
 
@@ -2396,7 +2398,7 @@ bool testAssignVar1D()
       h1->Fill(value, 1.0);
    }
 
-   TH1D* h2 = new TH1D("=1D-h2", "h2-Title", numberOfBins, v);
+   TH1D* h2 = new TH1D("=1D_h2", "h2-Title", numberOfBins, v);
    *h2 = *h1;
 
    bool ret = equals("Assign Oper VarH '='  1D", h1, h2, cmpOptStats);
@@ -2451,7 +2453,7 @@ bool testCopyConstructor1D()
 {
    // Tests the copy constructor for 1D Histograms
 
-   TH1D* h1 = new TH1D("cc1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("cc1D_h1", "h1-Title", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();
 
@@ -2474,7 +2476,7 @@ bool testCopyConstructorVar1D()
    Double_t v[numberOfBins+1];
    FillVariableRange(v);
 
-   TH1D* h1 = new TH1D("cc1D-h1", "h1-Title", numberOfBins, v);
+   TH1D* h1 = new TH1D("cc1D_h1", "h1-Title", numberOfBins, v);
 
    h1->Sumw2();
 
@@ -2535,7 +2537,7 @@ bool testClone1D()
 {
    // Tests the clone method for 1D Histograms
 
-   TH1D* h1 = new TH1D("cl1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("cl1D_h1", "h1-Title", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();
 
@@ -2558,7 +2560,7 @@ bool testCloneVar1D()
    Double_t v[numberOfBins+1];
    FillVariableRange(v);
 
-   TH1D* h1 = new TH1D("cl1D-h1", "h1-Title", numberOfBins, v);
+   TH1D* h1 = new TH1D("cl1D_h1", "h1-Title", numberOfBins, v);
 
    h1->Sumw2();
 
@@ -2619,7 +2621,7 @@ bool testAssign2D()
 {
    // Tests the operator=() method for 2D Histograms
 
-   TH2D* h1 = new TH2D("=2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("=2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -2631,7 +2633,7 @@ bool testAssign2D()
       h1->Fill(x, y, 1.0);
    }
 
-   TH2D* h2 = new TH2D("=2D-h2", "h2-Title",
+   TH2D* h2 = new TH2D("=2D_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    *h2 = *h1;
@@ -2671,7 +2673,7 @@ bool testCopyConstructor2D()
 {
    // Tests the copy constructor for 2D Histograms
 
-   TH2D* h1 = new TH2D("cc2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("cc2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -2716,7 +2718,7 @@ bool testClone2D()
 {
    // Tests the clone method for 2D Histograms
 
-   TH2D* h1 = new TH2D("cl2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("cl2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -2761,7 +2763,7 @@ bool testAssign3D()
 {
    // Tests the operator=() method for 3D Histograms
 
-   TH3D* h1 = new TH3D("=3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("=3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2775,7 +2777,7 @@ bool testAssign3D()
       h1->Fill(x, y, z, 1.0);
    }
 
-   TH3D* h2 = new TH3D("=3D-h2", "h2-Title",
+   TH3D* h2 = new TH3D("=3D_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2818,7 +2820,7 @@ bool testCopyConstructor3D()
 {
    // Tests the copy constructor for 3D Histograms
 
-   TH3D* h1 = new TH3D("cc3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("cc3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2867,7 +2869,7 @@ bool testClone3D()
 {
    // Tests the clone method for 3D Histograms
 
-   TH3D* h1 = new TH3D("cl3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("cl3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -2945,7 +2947,7 @@ bool testWriteRead1D()
 {
    // Tests the write and read methods for 1D Histograms
 
-   TH1D* h1 = new TH1D("wr1D-h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("wr1D_h1", "h1-Title", numberOfBins, minRange, maxRange);
 
    h1->Sumw2();
 
@@ -2959,7 +2961,7 @@ bool testWriteRead1D()
    f.Close();
 
    TFile f2("tmpHist.root");
-   TH1D* h2 = static_cast<TH1D*> ( f2.Get("wr1D-h1") );
+   TH1D* h2 = static_cast<TH1D*> ( f2.Get("wr1D_h1") );
 
    bool ret = equals("Read/Write Hist 1D", h1, h2, cmpOptStats);
    if (cleanHistos) delete h1;
@@ -2973,7 +2975,7 @@ bool testWriteReadVar1D()
    Double_t v[numberOfBins+1];
    FillVariableRange(v);
 
-   TH1D* h1 = new TH1D("wr1D-h1", "h1-Title", numberOfBins, v);
+   TH1D* h1 = new TH1D("wr1D_h1", "h1-Title", numberOfBins, v);
 
    h1->Sumw2();
 
@@ -2987,7 +2989,7 @@ bool testWriteReadVar1D()
    f.Close();
 
    TFile f2("tmpHist.root");
-   TH1D* h2 = static_cast<TH1D*> ( f2.Get("wr1D-h1") );
+   TH1D* h2 = static_cast<TH1D*> ( f2.Get("wr1D_h1") );
 
    bool ret = equals("Read/Write VarH 1D", h1, h2, cmpOptStats);
    if (cleanHistos) delete h1;
@@ -3049,7 +3051,7 @@ bool testWriteRead2D()
 {
    // Tests the write and read methods for 2D Histograms
 
-   TH2D* h1 = new TH2D("wr2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("wr2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -3066,7 +3068,7 @@ bool testWriteRead2D()
    f.Close();
 
    TFile f2("tmpHist.root");
-   TH2D* h2 = static_cast<TH2D*> ( f2.Get("wr2D-h1") );
+   TH2D* h2 = static_cast<TH2D*> ( f2.Get("wr2D_h1") );
 
    bool ret = equals("Read/Write Hist 2D", h1, h2, cmpOptStats);
    if (cleanHistos) delete h1;
@@ -3104,7 +3106,7 @@ bool testWriteRead3D()
 {
    // Tests the write and read methods for 3D Histograms
 
-   TH3D* h1 = new TH3D("wr3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("wr3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -3123,7 +3125,7 @@ bool testWriteRead3D()
    f.Close();
 
    TFile f2("tmpHist.root");
-   TH3D* h2 = static_cast<TH3D*> ( f2.Get("wr3D-h1") );
+   TH3D* h2 = static_cast<TH3D*> ( f2.Get("wr3D_h1") );
 
    bool ret = equals("Read/Write Hist 3D", h1, h2, cmpOptStats);
    if (cleanHistos) delete h1;
@@ -3350,16 +3352,16 @@ bool testMerge2D()
 {
    // Tests the merge method for 2D Histograms
 
-   TH2D* h1 = new TH2D("merge2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("merge2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("merge2D-h2", "h2-Title",
+   TH2D* h2 = new TH2D("merge2D_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h3 = new TH2D("merge2D-h3", "h3-Title",
+   TH2D* h3 = new TH2D("merge2D_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h4 = new TH2D("merge2D-h4", "h4-Title",
+   TH2D* h4 = new TH2D("merge2D_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -3457,19 +3459,19 @@ bool testMerge3D()
 {
    // Tests the merge method for 3D Histograms
 
-   TH3D* h1 = new TH3D("merge3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("merge3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("merge3D-h2", "h2-Title",
+   TH3D* h2 = new TH3D("merge3D_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("merge3D-h3", "h3-Title",
+   TH3D* h3 = new TH3D("merge3D_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h4 = new TH3D("merge3D-h4", "h4-Title",
+   TH3D* h4 = new TH3D("merge3D_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -3764,16 +3766,16 @@ bool testMerge2DLabelSame()
    // histogram with labels - just merges according to the x-values
    // This test is basically useless
 
-   TH2D* h1 = new TH2D("merge2DLabelSame-h1", "h1-Title",
+   TH2D* h1 = new TH2D("merge2DLabelSame_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("merge2DLabelSame-h2", "h2-Title",
+   TH2D* h2 = new TH2D("merge2DLabelSame_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h3 = new TH2D("merge2DLabelSame-h3", "h3-Title",
+   TH2D* h3 = new TH2D("merge2DLabelSame_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h4 = new TH2D("merge2DLabelSame-h4", "h4-Title",
+   TH2D* h4 = new TH2D("merge2DLabelSame_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -3820,19 +3822,19 @@ bool testMerge3DLabelSame()
 {
    // Tests the merge with some equal labels method for 3D Histograms
 
-   TH3D* h1 = new TH3D("merge3DLabelSame-h1", "h1-Title",
+   TH3D* h1 = new TH3D("merge3DLabelSame_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("merge3DLabelSame-h2", "h2-Title",
+   TH3D* h2 = new TH3D("merge3DLabelSame_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("merge3DLabelSame-h3", "h3-Title",
+   TH3D* h3 = new TH3D("merge3DLabelSame_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h4 = new TH3D("merge3DLabelSame-h4", "h4-Title",
+   TH3D* h4 = new TH3D("merge3DLabelSame_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -4065,10 +4067,10 @@ bool testMerge1DLabelDiff()
 {
    // Tests the merge with some different labels  for 1D Histograms
 
-   TH1D* h1 = new TH1D("merge1DLabelDiff-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("merge1DLabelDiff-h2", "h2-Title", numberOfBins, minRange, maxRange);
-   TH1D* h3 = new TH1D("merge1DLabelDiff-h3", "h3-Title", numberOfBins, minRange, maxRange);
-   TH1D* h4 = new TH1D("merge1DLabelDiff-h4", "h4-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("merge1DLabelDiff_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("merge1DLabelDiff_h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("merge1DLabelDiff_h3", "h3-Title", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("merge1DLabelDiff_h4", "h4-Title", numberOfBins, minRange, maxRange);
 
    // This test fails, as expected! That is why it is not run in the tests suite.
    const char labels[10][5] = {"aaa","bbb","ccc","ddd","eee","fff","ggg","hhh","iii","lll"};
@@ -4116,6 +4118,8 @@ bool testMerge1DLabelDiff()
    list->Add(h2);
    list->Add(h3);
 
+   if (!cleanHistos) h1->Clone("merge1DLabelDiff_h0");
+
    h1->Merge(list);
 
    // need to order the histo to compare them
@@ -4137,16 +4141,16 @@ bool testMerge2DLabelDiff()
    // are different ones and still the tests passes! This is not
    // consistent with TH1::Merge()
 
-   TH2D* h1 = new TH2D("merge2DLabelDiff-h1", "h1-Title",
+   TH2D* h1 = new TH2D("merge2DLabelDiff_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("merge2DLabelDiff-h2", "h2-Title",
+   TH2D* h2 = new TH2D("merge2DLabelDiff_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h3 = new TH2D("merge2DLabelDiff-h3", "h3-Title",
+   TH2D* h3 = new TH2D("merge2DLabelDiff_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h4 = new TH2D("merge2DLabelDiff-h4", "h4-Title",
+   TH2D* h4 = new TH2D("merge2DLabelDiff_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -4197,19 +4201,19 @@ bool testMerge3DLabelDiff()
    // are different ones and still the tests passes! This is not
    // consistent with TH1::Merge()
 
-   TH3D* h1 = new TH3D("merge3DLabelDiff-h1", "h1-Title",
+   TH3D* h1 = new TH3D("merge3DLabelDiff_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("merge3DLabelDiff-h2", "h2-Title",
+   TH3D* h2 = new TH3D("merge3DLabelDiff_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("merge3DLabelDiff-h3", "h3-Title",
+   TH3D* h3 = new TH3D("merge3DLabelDiff_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h4 = new TH3D("merge3DLabelDiff-h4", "h4-Title",
+   TH3D* h4 = new TH3D("merge3DLabelDiff_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -4501,16 +4505,16 @@ bool testMerge2DLabelAll()
 {
    // Tests the merge method with fully equally labelled 2D Histograms
 
-   TH2D* h1 = new TH2D("merge2DLabelAll-h1", "h1-Title",
+   TH2D* h1 = new TH2D("merge2DLabelAll_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("merge2DLabelAll-h2", "h2-Title",
+   TH2D* h2 = new TH2D("merge2DLabelAll_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h3 = new TH2D("merge2DLabelAll-h3", "h3-Title",
+   TH2D* h3 = new TH2D("merge2DLabelAll_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h4 = new TH2D("merge2DLabelAll-h4", "h4-Title",
+   TH2D* h4 = new TH2D("merge2DLabelAll_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
    for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
@@ -4560,19 +4564,19 @@ bool testMerge3DLabelAll()
 {
    // Tests the merge method with fully equally labelled 3D Histograms
 
-   TH3D* h1 = new TH3D("merge3DLabelAll-h1", "h1-Title",
+   TH3D* h1 = new TH3D("merge3DLabelAll_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("merge3DLabelAll-h2", "h2-Title",
+   TH3D* h2 = new TH3D("merge3DLabelAll_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("merge3DLabelAll-h3", "h3-Title",
+   TH3D* h3 = new TH3D("merge3DLabelAll_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h4 = new TH3D("merge3DLabelAll-h4", "h4-Title",
+   TH3D* h4 = new TH3D("merge3DLabelAll_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -4812,10 +4816,10 @@ bool testMerge1DLabelAllDiff()
    //LM: Dec 2010 : rmeake this test as
    // a test of histogram with some different labels not all filled
 
-   TH1D* h1 = new TH1D("merge1DLabelAllDiff-h1", "h1-Title", numberOfBins, minRange, maxRange);
-   TH1D* h2 = new TH1D("merge1DLabelAllDiff-h2", "h2-Title", numberOfBins, minRange, maxRange);
-   TH1D* h3 = new TH1D("merge1DLabelAllDiff-h3", "h3-Title", numberOfBins, minRange, maxRange);
-   TH1D* h4 = new TH1D("merge1DLabelAllDiff-h4", "h4-Title", numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("merge1DLabelAllDiff_h1", "h1-Title", numberOfBins, minRange, maxRange);
+   TH1D* h2 = new TH1D("merge1DLabelAllDiff_h2", "h2-Title", numberOfBins, minRange, maxRange);
+   TH1D* h3 = new TH1D("merge1DLabelAllDiff_h3", "h3-Title", numberOfBins, minRange, maxRange);
+   TH1D* h4 = new TH1D("merge1DLabelAllDiff_h4", "h4-Title", numberOfBins, minRange, maxRange);
 
    Int_t ibin = r.Integer(numberOfBins)+1;
    h1->GetXaxis()->SetBinLabel(ibin,"aaa");
@@ -4873,17 +4877,17 @@ bool testMerge2DLabelAllDiff()
    // Note: in case of underflow/overflow in x axis not clear  how merge should proceed
    // when merging with labels underflow/overflow will not be considered
 
-   TH2D* h1 = new TH2D("merge2DLabelAllDiff-h1", "h1-Title",
+   TH2D* h1 = new TH2D("merge2DLabelAllDiff_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange);
-   TH2D* h2 = new TH2D("merge2DLabelAllDiff-h2", "h2-Title",
+   TH2D* h2 = new TH2D("merge2DLabelAllDiff_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange);
-   TH2D* h3 = new TH2D("merge2DLabelAllDiff-h3", "h3-Title",
+   TH2D* h3 = new TH2D("merge2DLabelAllDiff_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange);
-   TH2D* h4 = new TH2D("merge2DLabelAllDiff-h4", "h4-Title",
-                       numberOfBins, minRange, maxRange,
+   TH2D* h4 = new TH2D("merge2DLabelAllDiff_h4", "h4-Title",
+                       2*numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange);
 
 
@@ -4900,6 +4904,9 @@ bool testMerge2DLabelAllDiff()
       // use for h3 same label as for h2 to test the merging
       h3->GetXaxis()->SetBinLabel(i, name2.str().c_str());
       h3->GetYaxis()->SetBinLabel(i, name2.str().c_str());
+       // we set the bin labels also in h4
+      h4->GetXaxis()->SetBinLabel(i, name1.str().c_str());
+      h4->GetXaxis()->SetBinLabel(i+numberOfBins, name2.str().c_str());
    }
 
    // the x axis will be full labels while the y axis will be numeric
@@ -4955,19 +4962,19 @@ bool testMerge3DLabelAllDiff()
    // All label sizes are less than number of bins, therefore axis cannot be extended
    // and merge is done then numerically and not in label mode
 
-   TH3D* h1 = new TH3D("merge3DLabelAllDiff-h1", "h1-Title",
+   TH3D* h1 = new TH3D("merge3DLabelAllDiff_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("merge3DLabelAllDiff-h2", "h2-Title",
+   TH3D* h2 = new TH3D("merge3DLabelAllDiff_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h3 = new TH3D("merge3DLabelAllDiff-h3", "h3-Title",
+   TH3D* h3 = new TH3D("merge3DLabelAllDiff_h3", "h3-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h4 = new TH3D("merge3DLabelAllDiff-h4", "h4-Title",
+   TH3D* h4 = new TH3D("merge3DLabelAllDiff_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -5311,16 +5318,16 @@ bool testMerge2D_Diff(bool testEmpty = false)
    //LM. t.b.u.: for 1D can make h3 with 330 bins , while in 2D if I make h3 with 33 bins
    //  routine which check axis fails. Needs to be improved ???
 
-   TH2D *h1 = new TH2D("merge2DDiff-h1","h1-Title",
+   TH2D *h1 = new TH2D("merge2DDiff_h1","h1-Title",
                        11,-110,0,
                        11,-110,0);
-   TH2D *h2 = new TH2D("merge2DDiff-h2","h2-Title",
+   TH2D *h2 = new TH2D("merge2DDiff_h2","h2-Title",
                        22,0,110,
                        22,0,110);
-   TH2D *h3 = new TH2D("merge2DDiff-h3","h3-Title",
+   TH2D *h3 = new TH2D("merge2DDiff_h3","h3-Title",
                        44,-55,55,
                        44,-55,55);
-   TH2D *h4 = new TH2D("merge2DDiff-h4","h4-Title",
+   TH2D *h4 = new TH2D("merge2DDiff_h4","h4-Title",
                        22,-110,110,
                        22,-110,110);
 
@@ -5375,19 +5382,19 @@ bool testMerge3D_Diff(bool testEmpty = false)
    // Tests the merge method with different binned 3D Histograms
 
 
-   TH3D *h1 = new TH3D("merge3DDiff-h1","h1-Title",
+   TH3D *h1 = new TH3D("merge3DDiff_h1","h1-Title",
                        11,-110,0,
                        11,-110,0,
                        11,-110,0);
-   TH3D *h2 = new TH3D("merge3DDiff-h2","h2-Title",
+   TH3D *h2 = new TH3D("merge3DDiff_h2","h2-Title",
                        22,0,110,
                        22,0,110,
                        22,0,110);
-   TH3D *h3 = new TH3D("merge3DDiff-h3","h3-Title",
+   TH3D *h3 = new TH3D("merge3DDiff_h3","h3-Title",
                        44,-55,55,
                        44,-55,55,
                        44,-55,55);
-   TH3D *h4 = new TH3D("merge3DDiff-h4","h4-Title",
+   TH3D *h4 = new TH3D("merge3DDiff_h4","h4-Title",
                        22,-110,110,
                        22,-110,110,
                        22,-110,110);
@@ -5657,13 +5664,13 @@ bool testMerge2DExtend(UInt_t extendType = TH1::kAllAxes)
    // Tests the merge method for diferent 1D Histograms
    // when axis can be extended (e.g. for time histograms)
 
-   TH2D* h1 = new TH2D("merge2D-h1", "h1-Title",
+   TH2D* h1 = new TH2D("merge2D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h2 = new TH2D("merge2D-h2", "h2-Title",
+   TH2D* h2 = new TH2D("merge2D_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH2D* h4 = new TH2D("merge2D-h4", "h4-Title",
+   TH2D* h4 = new TH2D("merge2D_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
 
@@ -5714,15 +5721,15 @@ bool testMerge3DExtend(UInt_t extendType = TH1::kAllAxes)
    // Tests the merge method for diferent 1D Histograms
    // when axis can be extended (e.g. for time histograms)
 
-   TH3D* h1 = new TH3D("merge3D-h1", "h1-Title",
+   TH3D* h1 = new TH3D("merge3D_h1", "h1-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h2 = new TH3D("merge3D-h2", "h2-Title",
+   TH3D* h2 = new TH3D("merge3D_h2", "h2-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
-   TH3D* h4 = new TH3D("merge3D-h4", "h4-Title",
+   TH3D* h4 = new TH3D("merge3D_h4", "h4-Title",
                        numberOfBins, minRange, maxRange,
                        numberOfBins + 1, minRange, maxRange,
                        numberOfBins + 2, minRange, maxRange);
@@ -5883,19 +5890,27 @@ bool testLabel()
 {
    // Tests labelling a 1D Histogram
 
-   TH1D* h1 = new TH1D("lD1-h1", "h1-Title", 2*numberOfBins, minRange, maxRange);
+   TH1D* h1 = new TH1D("lD1_h1", "h1-Title", 2*numberOfBins, minRange, maxRange);
    // build histo with extra  labels to tets the deflate option
    int extraBins = 20;
-   TH1D* h2 = new TH1D("lD1-h2", "h2-Title", 2*numberOfBins+20, minRange, maxRange + extraBins*h1->GetXaxis()->GetBinWidth(1));
+   TH1D* h2 = new TH1D("lD1_h2", "h2-Title", 2*numberOfBins+20, minRange, maxRange + extraBins*h1->GetXaxis()->GetBinWidth(1));
 
 
    // set labels
-   std::vector<std::string> vLabels;
-   for ( Int_t bin = 1; bin <= h1->GetNbinsX() ; ++bin ) {
+   std::vector<std::string> vLabels(h1->GetNbinsX());
+   std::vector<int> bins(h1->GetNbinsX());
+   for ( Int_t i = 0; i < h1->GetNbinsX() ; ++i ) {
+      Int_t bin = i+1;
       ostringstream label;
-      label << bin;
-      vLabels.push_back(label.str());
-      h2->GetXaxis()->SetBinLabel(bin, label.str().c_str());
+      char letter = (char) ((int) 'a' + i );
+      label << letter;
+      vLabels[i] = label.str();
+      bins[i] = bin;
+   }
+   // set bin label in random order in bins to test ordering when labels are filled randomly
+   std::shuffle(bins.begin(), bins.end(), std::default_random_engine{});
+   for (size_t i = 0; i < bins.size(); ++i ) {
+      h2->GetXaxis()->SetBinLabel(bins[i], vLabels[i].c_str());
    }
    // sort labels in alphabetic order
    std::sort(vLabels.begin(), vLabels.end() );
@@ -5909,10 +5924,13 @@ bool testLabel()
       h2->Fill(vLabels[bin-1].c_str(), 1.0);
    }
 
+   // test ordering label in content ascending order
+   //h2->LabelsOption("<","x");
+   // test ordering label alphabetically
    h2->LabelsOption("a");
    h2->LabelsDeflate();
 
-   bool status = equals("Fill(char*)", h1, h2, cmpOptStats, 1E-13);
+   bool status = equals("testLabel", h1, h2, cmpOptStats, 1E-13);
    if (cleanHistos) delete h1;
    return status;
 }
@@ -5922,21 +5940,32 @@ bool testLabel2DX()
 {
    // Tests labelling a 1D Histogram
 
-   TH2D* h1 = new TH2D("lD2-h1", "h1-Title", 2*numberOfBins, minRange, maxRange, numberOfBins, minRange, maxRange);
+   TH2D* h1 = new TH2D("lD2_h1", "h1-Title", 2*numberOfBins, minRange, maxRange, numberOfBins, minRange, maxRange);
    // build histo with extra  labels to tets the deflate option
-   TH2D* h2 = new TH2D("lD2-h2", "h2-Title", 2*numberOfBins+20, minRange, maxRange + 20*h1->GetXaxis()->GetBinWidth(1), numberOfBins, minRange, maxRange);
+   TH2D* h2 = new TH2D("lD2_h2", "h2-Title", 2*numberOfBins+20, minRange, maxRange + 20*h1->GetXaxis()->GetBinWidth(1), numberOfBins, minRange, maxRange);
 
    // set labels
-   std::vector<std::string> vLabels;
-   for ( Int_t bin = 1; bin <= h1->GetNbinsX() ; ++bin ) {
+   std::vector<std::string> vLabels(h1->GetNbinsX());
+   std::vector<int> bins(h1->GetNbinsX());
+   for ( Int_t i = 0; i < h1->GetNbinsX() ; ++i ) {
+      Int_t bin = i+1;
       ostringstream label;
-      label << bin;
-      vLabels.push_back(label.str());
-      h2->GetXaxis()->SetBinLabel(bin, label.str().c_str());
+      char letter = (char) ((int) 'a' + i );
+      label << letter;
+      vLabels[i] = label.str();
+      bins[i] = bin;
+   }
+   // set bin label in random order in bins to test ordering when labels are filled randomly
+   std::shuffle(bins.begin(), bins.end(), std::default_random_engine{});
+   for (size_t i = 0; i < bins.size(); ++i ) {
+      h2->GetXaxis()->SetBinLabel(bins[i], vLabels[i].c_str());
    }
    // sort labels in alphabetic order
    std::sort(vLabels.begin(), vLabels.end() );
 
+   // fill h1 with numbers and h2 using labels
+   // since labels are ordered alphabetically
+   // for filling bin i-th of h1 same bin i-th will be filled of h2
    for ( Int_t e = 0; e < nEvents; ++e ) {
       Double_t xvalue = r.Uniform(minRange, maxRange);
       Double_t yvalue = r.Uniform(minRange, maxRange);
@@ -5947,11 +5976,14 @@ bool testLabel2DX()
       h2->Fill( vLabels[binx-1].c_str(), h1->GetYaxis()->GetBinCenter(biny), 1.0);
    }
 
-   h2->LabelsOption("a");
+   // test ordering label in content descending order
+   //h2->LabelsOption(">","x");
+   // test ordering label alphabetically
+   h2->LabelsOption("a","x");
 
    h2->LabelsDeflate();
 
-   bool status = equals("Fill(char*)", h1, h2, cmpOptStats, 1E-13);
+   bool status = equals("testLabel2DX", h1, h2, cmpOptStats, 1E-13);
    if (cleanHistos) delete h1;
    return status;
 }
@@ -5960,9 +5992,9 @@ bool testLabel2DY()
 {
    // Tests labelling a 1D Histogram
 
-   TH2D* h1 = new TH2D("lD2-h1", "h1-Title", numberOfBins, minRange, maxRange, 2*numberOfBins, minRange, maxRange);
+   TH2D* h1 = new TH2D("lD2_h1", "h1-Title", numberOfBins, minRange, maxRange, 2*numberOfBins, minRange, maxRange);
    // build histo with extra  labels to tets the deflate option
-   TH2D* h2 = new TH2D("lD2-h2", "h2-Title", numberOfBins, minRange, maxRange, 2*numberOfBins+20, minRange, maxRange + 20*h1->GetYaxis()->GetBinWidth(1));
+   TH2D* h2 = new TH2D("lD2_h2", "h2-Title", numberOfBins, minRange, maxRange, 2*numberOfBins+20, minRange, maxRange + 20*h1->GetYaxis()->GetBinWidth(1));
 
    // set labels
    std::vector<std::string> vLabels;
@@ -5985,11 +6017,11 @@ bool testLabel2DY()
       h2->Fill(  h1->GetXaxis()->GetBinCenter(binx), vLabels[biny-1].c_str(), 1.0);
    }
 
-   h2->GetYaxis()->LabelsOption("a");
+   h2->LabelsOption("a","y");
 
    h2->LabelsDeflate("Y");
 
-   bool status = equals("Fill(char*)", h1, h2, cmpOptStats, 1E-13);
+   bool status = equals("testLabel2DY", h1, h2, cmpOptStats, 1E-13);
    if (cleanHistos) delete h1;
    return status;
 }
@@ -7752,6 +7784,10 @@ bool testRefRead1D()
       h1->Write();
    } else {
       h1 = static_cast<TH1D*> ( refFile->Get("rr1D-h1") );
+      if (!h1) {
+          Error("testRefRead1D","Error reading histogram rr1D-h1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       TH1D* h2 = new TH1D("rr1D-h2", "h2-Title", numberOfBins, minRange, maxRange);
       h2->Sumw2();
 
@@ -7786,6 +7822,10 @@ bool testRefReadProf1D()
    } else {
       TH1::SetDefaultSumw2(false);
       p1 = static_cast<TProfile*> ( refFile->Get("rr1D-p1") );
+      if (!p1) {
+          Error("testRefReadProf1D","Error reading profile rr1D-p1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       TProfile* p2 = new TProfile("rr1D-p2", "p2-Title", numberOfBins, minRange, maxRange);
 //      p2->Sumw2();
 
@@ -7823,6 +7863,10 @@ bool testRefRead2D()
       h1->Write();
    } else {
       h1 = static_cast<TH2D*> ( refFile->Get("rr2D-h1") );
+      if (!h1) {
+          Error("testRefRead2D","Error reading histogram rr2D-h1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       TH2D* h2 = new TH2D("rr2D-h2", "h2-Title",
                           numberOfBins, minRange, maxRange,
                           numberOfBins, minRange, maxRange);
@@ -7860,6 +7904,10 @@ bool testRefReadProf2D()
       p1->Write();
    } else {
       p1 = static_cast<TProfile2D*> ( refFile->Get("rr2D-p1") );
+      if (!p1) {
+          Error("testRefReadProf2D","Error reading profile rr2D-p1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       TProfile2D* p2 = new TProfile2D("rr2D-p2", "p2-Title",
                                       numberOfBins, minRange, maxRange,
                                       numberOfBins, minRange, maxRange);
@@ -7899,6 +7947,10 @@ bool testRefRead3D()
       h1->Write();
    } else {
       h1 = static_cast<TH3D*> ( refFile->Get("rr3D-h1") );
+      if (!h1) {
+          Error("testRefRead3D","Error reading histogram rr3D-h1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       TH3D* h2 = new TH3D("rr3D-h2", "h2-Title",
                           numberOfBins, minRange, maxRange,
                           numberOfBins, minRange, maxRange,
@@ -7940,6 +7992,10 @@ bool testRefReadProf3D()
       p1->Write();
    } else {
       p1 = static_cast<TProfile3D*> ( refFile->Get("rr3D-p1") );
+      if (!p1) {
+          Error("testRefReadProf3D","Error reading profile rr3D-p1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       TProfile3D* p2 = new TProfile3D("rr3D-p2", "p2-Title",
                           numberOfBins, minRange, maxRange,
                           numberOfBins, minRange, maxRange,
@@ -7986,6 +8042,10 @@ bool testRefReadSparse()
       s1->Write();
    } else {
       s1 = static_cast<THnSparseD*> ( refFile->Get("rr-s1") );
+      if (!s1) {
+          Error("testRefReadSparse","Error reading THnSparse rr-s1 from file");
+          return kTRUE;  // true indicates a failure
+      }
       THnSparseD* s2 = new THnSparseD("rr-s1", "s1-Title", 3, bsize, xmin, xmax);
       s2->Sumw2();
 
@@ -10328,7 +10388,7 @@ int stressHistogram(int testNumber = 0)
                                           "Reference File Read for Histograms and Profiles..................",
                                           refReadTestPointer };
 
-
+   // test24 - compare with a reference old file
    testCounter++;
    if (runAll || testNumber == testCounter) {
       if (refFileOption == refFileWrite) {
@@ -10343,6 +10403,10 @@ int stressHistogram(int testNumber = 0)
 
       if (refFile != 0) {
          r.SetSeed(8652);
+         if (defaultEqualOptions == cmpOptDebug) {
+            std::cout << "content of file " << refFile->GetName() << std::endl;
+            refFile->ls();
+         }
          status = 0;
          for (unsigned int j = 0; j < refReadTestSuite.nTests; ++j) {
             status += refReadTestSuite.tests[j]();
@@ -10606,6 +10670,9 @@ int equals(const char* msg, TH2D* h1, TH2D* h2, int options, double ERRORLIMIT)
            << (h1 == h2 ) << " " << differents << std::endl;
    }
 
+   bool labelXaxis = (h1->GetXaxis()->GetLabels() && h1->GetXaxis()->CanExtend());
+   bool labelYaxis = (h1->GetYaxis()->GetLabels() && h1->GetYaxis()->CanExtend());
+
    for ( int i = 0; i <= h1->GetNbinsX() + 1; ++i )
       for ( int j = 0; j <= h1->GetNbinsY() + 1; ++j )
       {
@@ -10614,9 +10681,15 @@ int equals(const char* msg, TH2D* h1, TH2D* h2, int options, double ERRORLIMIT)
 
          if (debug)
          {
-            std::cout << equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT) << " "
-                 << equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT) << " "
-                 << "[" << x << "," << y << "]: "
+            if (!labelXaxis)
+               std::cout << equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT) << " ";
+            else
+               std::cout << equals(h1->GetXaxis()->GetBinLabel(i), h2->GetXaxis()->GetBinLabel(i) ) << " ";
+            if (!labelYaxis)
+               std::cout << equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT) << " ";
+            else
+               std::cout << equals(h1->GetYaxis()->GetBinLabel(j), h2->GetYaxis()->GetBinLabel(j) ) << " ";
+            std::cout  << "[" << i << " : " << x << ", " << j << " : " << y << "]: "
                  << h1->GetBinContent(i,j) << " +/- " << h1->GetBinError(i,j) << " | "
                  << h2->GetBinContent(i,j) << " +/- " << h2->GetBinError(i,j)
                  << " | " << equals(h1->GetBinContent(i,j), h2->GetBinContent(i,j), ERRORLIMIT)
@@ -10625,8 +10698,15 @@ int equals(const char* msg, TH2D* h1, TH2D* h2, int options, double ERRORLIMIT)
                  << " "   << (fabs(h1->GetBinContent(i,j) - h2->GetBinContent(i,j)))
                  << std::endl;
          }
-         differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
-         differents += (bool) equals(y, h2->GetYaxis()->GetBinCenter(j), ERRORLIMIT);
+         if (labelXaxis)
+            differents += equals(h1->GetXaxis()->GetBinLabel(i), h2->GetXaxis()->GetBinLabel(i) );
+         else
+            differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
+         if (labelYaxis)
+            differents += equals(h1->GetYaxis()->GetBinLabel(j), h2->GetYaxis()->GetBinLabel(j) );
+         else
+            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);
@@ -10657,28 +10737,30 @@ int equals(const char* msg, TH1D* h1, TH1D* h2, int options, double ERRORLIMIT)
            << (h1 == h2 ) << " " << differents << std::endl;
    }
 
-   // check axis
+   bool labelAxis = (h1->GetXaxis()->GetLabels() && h1->GetXaxis()->CanExtend());
 
    differents += (bool) equals(h1->GetXaxis()->GetNbins() , h2->GetXaxis()->GetNbins() );
    if (debug) {
       cout << "Nbins  = " << h1->GetXaxis()->GetNbins() << " |  " <<  h2->GetXaxis()->GetNbins() << " | " << differents << std::endl;
    }
 
-   differents += (bool) equals(h1->GetXaxis()->GetXmin() , h2->GetXaxis()->GetXmin() );
-   if (debug) {
-      cout << "Xmin   = "  << h1->GetXaxis()->GetXmin() << " |  " <<  h2->GetXaxis()->GetXmin() << " | " << differents << std::endl;
-   }
+   if (!labelAxis) {
+      differents += (bool) equals(h1->GetXaxis()->GetXmin() , h2->GetXaxis()->GetXmin() );
+      if (debug) {
+         cout << "Xmin   = "  << h1->GetXaxis()->GetXmin() << " |  " <<  h2->GetXaxis()->GetXmin() << " | " << differents << std::endl;
+      }
 
-   differents += (bool) equals(h1->GetXaxis()->GetXmax() , h2->GetXaxis()->GetXmax() );
-   if (debug) {
-      cout << "Xmax   = "  << h1->GetXaxis()->GetXmax() << " |  " <<  h2->GetXaxis()->GetXmax() << endl;
+      differents += (bool) equals(h1->GetXaxis()->GetXmax() , h2->GetXaxis()->GetXmax() );
+      if (debug) {
+         cout << "Xmax   = "  << h1->GetXaxis()->GetXmax() << " |  " <<  h2->GetXaxis()->GetXmax() << endl;
+      }
    }
 
    for ( int i = 0; i <= h1->GetNbinsX() + 1; ++i )
    {
       Double_t x = h1->GetXaxis()->GetBinCenter(i);
 
-      differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
+      if (!labelAxis) differents += (bool) equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT);
       differents += (bool) equals(h1->GetBinContent(i), h2->GetBinContent(i), ERRORLIMIT);
 
       if ( compareError )
@@ -10687,7 +10769,7 @@ int equals(const char* msg, TH1D* h1, TH1D* h2, int options, double ERRORLIMIT)
       if ( debug )
       {
          std::cout << equals(x, h2->GetXaxis()->GetBinCenter(i), ERRORLIMIT)
-              << " [" << x << "]: "
+              << " [" << i << " : " << x << "]: "
               << h1->GetBinContent(i) << " +/- " << h1->GetBinError(i) << " | "
               << h2->GetBinContent(i) << " +/- " << h2->GetBinError(i)
               << " | " << equals(h1->GetBinContent(i), h2->GetBinContent(i), ERRORLIMIT)
@@ -10716,6 +10798,13 @@ int equals(Double_t n1, Double_t n2, double ERRORLIMIT)
       return fabs(n2) > ERRORLIMIT;
 }
 
+int equals(const char * s1, const char * s2)
+{
+   std::string name1(s1);
+   std::string name2(s2);
+   return name1 != name2;
+}
+
 int compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT)
 {
    int differents = 0;
@@ -10825,7 +10914,7 @@ int main(int argc, char** argv)
          defaultEqualOptions = cmpOptDebug;
       } else if (arg == "-fast") {
          cout << "stressHistogram: running in fast mode " << endl;
-         nEvents = 10;
+         nEvents = 20;
       } else if (arg == "-n") {
          cout << "stressHistogram: running single test" << endl;
          testNumber = atoi(argv[++i]);