stressHistogram.cxx

// @(#)root/test:$name:  $:$id: stressHistogram.cxx,v 1.15 2002/10/25 10:47:51 rdm exp $
// Authors: David Gonzalez Maline November 2008

//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*//
//                                                                               //
//                                                                               //
// Here there is a set of tests for the histogram classes (including             //
// histograms and profiles). The methods tested work on:                         //
//                                                                               //
// 1. Projection testing (with and without weights)                              //
// 2. Rebinning                                                                  //
// 3. Addition, multiplication an division operations.                           //
// 4. Building and copying instances.                                            //
// 5. I/O functionality (including reference with older versions).               //
// 6. Labeling.                                                                  //
// 7. Interpolation                                                              //
//                                                                               //
// To see the tests individually, at the bottom of the file the tests            //
// are exectued using the structure TTestSuite, that defines the                 //
// subset, the number of routines to be tested as well as the pointes            //
// for these. Every tests is mean to be simple enough to be understood           //
// without much comments.                                                        //
//                                                                               //
// Finally, for debugging reasons, the struct compareOptions can be              //
// used to define the level of output of the tests, beging set                   //
// generally for the whole suit in defaultEqualOptions.                          //
// >> stressHistogram 1      : to print result for all tests                     //
// >> stressHistogram 2      : ro print each comparison, done for each bin       //
//                                                                               //
// An example of output when all the tests run OK is shown below:                //
// ****************************************************************************
// *  Starting  stress  H I S T O G R A M                                     *
// ****************************************************************************
// Test  1: Testing Histogram Projections without weights....................OK
// Test  2: Testing Profile Projections without weights......................OK
// Test  3: Testing Histogram Projections with weights.......................OK
// Test  4: Testing Profile   Projections with weights.......................OK
// Test  5: Projection with Range for Histograms and Profiles................OK
// Test  6: Histogram Rebinning..............................................OK
// Test  7: Add tests for 1D, 2D and 3D Histograms and Profiles..............OK
// Test  8: Multiply tests for 1D, 2D and 3D Histograms......................OK
// Test  9: Divide tests for 1D, 2D and 3D Histograms........................OK
// Test 10: Copy tests for 1D, 2D and 3D Histograms and Profiles.............OK
// Test 11: Read/Write tests for 1D, 2D and 3D Histograms and Profiles.......OK
// Test 12: Merge tests for 1D, 2D and 3D Histograms and Profiles............OK
// Test 13: Label tests for 1D and 2D Histograms ............................OK
// Test 14: Interpolation tests for Histograms...............................OK
// Test 15: Scale tests for Profiles.........................................OK
// Test 16: Integral tests for Histograms....................................OK
// Test 17: Buffer tests for Histograms......................................OK
// Test 18: Extend axis tests for Histograms.................................OK
// Test 19: TH1-THn[Sparse] Conversion tests.................................OK
// Test 20: FillData tests for Histograms and Sparses........................OK
// Test 21: Reference File Read for Histograms and Profiles..................OK
// ****************************************************************************
// stressHistogram: Real Time =  86.22 seconds Cpu Time =  85.64 seconds
//  ROOTMARKS = 1292.62 ROOT version: 6.05/01      remotes/origin/master@v6-05-01-336-g5c3d5ff
// ****************************************************************************
//                                                                               //
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*//


#include <sstream>
#include <cmath>

#include "TH2.h"
#include "TH3.h"
#include "TH2.h"
#include "THn.h"
#include "THnSparse.h"

#include "TProfile.h"
#include "TProfile2D.h"
#include "TProfile3D.h"

#include "TF1.h"
#include "TF2.h"
#include "TF3.h"

#include "Fit/SparseData.h"
#include "HFitInterface.h"
#include "TFitResult.h"

#include "Math/IntegratorOptions.h"

#include "TApplication.h"
#include "TBenchmark.h"
#include "Riostream.h"
#include "TMath.h"
#include "TRandom2.h"
#include "TFile.h"
#include "TClass.h"
#include "THashList.h"

#include "TROOT.h"
#include <algorithm>
#include <cassert>

using namespace std;

const unsigned int __DRAW__ = 0;

Double_t minRange = 1;
Double_t maxRange = 5;

const Double_t minRebin = 3;
const Double_t maxRebin = 7;

int nEvents = 1000;
const int numberOfBins = 10;

enum compareOptions {
   cmpOptNone=0,
   cmpOptPrint=1,
   cmpOptDebug=2,
   cmpOptNoError=4,
   cmpOptStats=8
};

int defaultEqualOptions = 0; //cmpOptPrint;
//int defaultEqualOptions = cmpOptDebug;

constexpr Bool_t cleanHistos = kTRUE;   // delete histogram after testing (swicth off in case of debugging)

const double defaultErrorLimit = 1.E-10;

enum RefFileEnum {
   refFileRead = 1,
   refFileWrite = 2
};

const int refFileOption = 1;
TFile * refFile = 0;
const char* refFileName = "http://root.cern.ch/files/stressHistogram.5.18.00.root";

TRandom2 r;
// set to zero if want to run different numbers every time
const int initialSeed = 0;



typedef bool ( * pointer2Test) ();

struct TTestSuite {
   unsigned int nTests;
   char suiteName[75];
   pointer2Test* tests;
};

// Methods for histogram comparisions (later implemented)
void printResult(const char* msg, bool status);
void FillVariableRange(Double_t v[numberOfBins+1]);
void FillHistograms(TH1D* h1, TH1D* h2, Double_t c1 = 1.0, Double_t c2 = 1.0);
void FillProfiles(TProfile* p1, TProfile* p2, Double_t c1 = 1.0, Double_t c2 = 1.0);
int equals(const char* msg, TH1D* h1, TH1D* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
int equals(const char* msg, TH2D* h1, TH2D* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
int equals(const char* msg, TH3D* h1, TH3D* h2, int options = 0, double ERRORLIMIT = defaultErrorLimit);
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 compareStatistics( TH1* h1, TH1* h2, bool debug, double ERRORLIMIT = defaultErrorLimit);
std::ostream& operator<<(std::ostream& out, TH1D* h);
// old stresHistOpts.cxx file

bool testAdd1()
{
   // Tests the first Add method for 1D Histograms

   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);

   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);
   h4->Add(h1, h2, c1, c2);

   bool ret = equals("Add1D1", h3, h4, cmpOptStats, 1E-13);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testAddProfile1()
{
   // Tests the first Add method for 1D Profiles

   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()
{
   // Tests the second Add method for 1D Histograms

   Double_t c2 = r.Rndm();

   TH1D* h5 = new TH1D("t1D2-h5", "h5=   h6+c2*h7", numberOfBins, minRange, maxRange);
   TH1D* h6 = new TH1D("t1D2-h6", "h6-Title", numberOfBins, minRange, maxRange);
   TH1D* h7 = new TH1D("t1D2-h7", "h7-Title", numberOfBins, minRange, maxRange);

   h5->Sumw2();h6->Sumw2();h7->Sumw2();

   FillHistograms(h6, h5, 1.0, 1.0);
   FillHistograms(h7, h5, 1.0, c2);

   h6->Add(h7, c2);

   bool ret = equals("Add1D2", h5, h6, cmpOptStats, 1E-13);
   if (cleanHistos) delete h5;
   if (cleanHistos) delete h7;
   return ret;
}

bool testAddProfile2()
{
   // Tests the second Add method for 1D Profiles

   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 testAdd3()
{
   // Tests the first add method to do scalation of 1D Histograms

   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);

   h1->Sumw2();h2->Sumw2();

   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value,  1.0);
      h2->Fill(value, c1 / h1->GetBinWidth( h1->FindBin(value) ) );
   }


   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
   h2->ResetStats();

   bool ret = equals("Add1D3", h2, h3, cmpOptStats, 1E-13);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   return ret;
}

bool testAddVar1()
{
   // Tests the second Add method for 1D Histograms with variable bin size

   Double_t v[numberOfBins+1];
   FillVariableRange(v);

   Double_t c1 = r.Rndm();
   Double_t c2 = r.Rndm();

   TH1D* h1 = new TH1D("h1", "h1-Title", numberOfBins, v);
   TH1D* h2 = new TH1D("h2", "h2-Title", numberOfBins, v);
   TH1D* h3 = new TH1D("h3", "h3=c1*h1+c2*h2", numberOfBins, v);

   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, v);
   h4->Add(h1, h2, c1, c2);

   bool ret = equals("AddVar1D1", h3, h4, cmpOptStats, 1E-13);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testAddVarProf1()
{

   // Tests the first Add method for 1D Profiles with variable bin size

   Double_t v[numberOfBins+1];
   FillVariableRange(v);

   Double_t c1 = r.Rndm();
   Double_t c2 = r.Rndm();

   TProfile* p1 = new TProfile("t1D1-p1", "p1-Title", numberOfBins, v);
   TProfile* p2 = new TProfile("t1D1-p2", "p2-Title", numberOfBins, v);
   TProfile* p3 = new TProfile("t1D1-p3", "p3=c1*p1+c2*p2", numberOfBins, v);

   FillProfiles(p1, p3, 1.0, c1);
   FillProfiles(p2, p3, 1.0, c2);

   TProfile* p4 = new TProfile("t1D1-p4", "p4=c1*p1+p2*c2", numberOfBins, v);
   p4->Add(p1, p2, c1, c2);

   bool ret = equals("AddVar1DProf1", p3, p4, cmpOptStats, 1E-13);
   delete p1;
   delete p2;
   delete p3;

   return ret;
}

bool testAddVar2()
{
   // Tests the second Add method for 1D Histograms with variable bin size

   Double_t v[numberOfBins+1];
   FillVariableRange(v);

   Double_t c2 = r.Rndm();

   TH1D* h5 = new TH1D("t1D2-h5", "h5=   h6+c2*h7", numberOfBins, v);
   TH1D* h6 = new TH1D("t1D2-h6", "h6-Title", numberOfBins, v);
   TH1D* h7 = new TH1D("t1D2-h7", "h7-Title", numberOfBins, v);

   h5->Sumw2();h6->Sumw2();h7->Sumw2();

   FillHistograms(h6, h5, 1.0, 1.0);
   FillHistograms(h7, h5, 1.0, c2);

   h6->Add(h7, c2);

   bool ret = equals("AddVar1D2", h5, h6, cmpOptStats, 1E-13);
   if (cleanHistos) delete h5;
   if (cleanHistos) delete h7;
   return ret;
}

bool testAddVarProf2()
{
   // Tests the second Add method for 1D Profiles with variable bin size

   Double_t v[numberOfBins+1];
   FillVariableRange(v);

   Double_t c2 = r.Rndm();

   TProfile* p5 = new TProfile("t1D2-p5", "p5=   p6+c2*p7", numberOfBins, v);
   TProfile* p6 = new TProfile("t1D2-p6", "p6-Title", numberOfBins, v);
   TProfile* p7 = new TProfile("t1D2-p7", "p7-Title", numberOfBins, v);

   p5->Sumw2();p6->Sumw2();p7->Sumw2();

   FillProfiles(p6, p5, 1.0, 1.0);
   FillProfiles(p7, p5, 1.0, c2);

   p6->Add(p7, c2);

   bool ret = equals("AddVar1D2", p5, p6, cmpOptStats, 1E-13);
   delete p5;
   delete p7;
   return ret;
}

bool testAddVar3()
{
   // Tests the first add method to do scale of 1D Histograms with variable bin width

   Double_t v[numberOfBins+1];
   FillVariableRange(v);

   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);

   h1->Sumw2();h2->Sumw2();

   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value,  1.0);
      h2->Fill(value, c1 / h1->GetBinWidth( h1->FindBin(value) ) );
   }

   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
   h2->ResetStats();

   bool ret = equals("Add1D3", h2, h3, cmpOptStats, 1E-13);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   return ret;
}


bool testAdd2D3()
{
   // Tests the first add method to do scale of 2D Histograms

   Double_t c1 = r.Rndm();

   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",
                       numberOfBins, minRange, maxRange,
                       numberOfBins+2, minRange, maxRange);

   h1->Sumw2();h2->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);
      Int_t binx = h1->GetXaxis()->FindBin(x);
      Int_t biny = h1->GetYaxis()->FindBin(y);
      Double_t area = h1->GetXaxis()->GetBinWidth( binx ) * h1->GetYaxis()->GetBinWidth( biny );
      h2->Fill(x, y, c1 / area);
   }

   TH2D* h3 = new TH2D("t1D1-h3", "h3=c1*h1",
                       numberOfBins, minRange, maxRange,
                       numberOfBins+2, 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
   h2->ResetStats();

   bool ret = equals("Add1D2", h2, h3, cmpOptStats, 1E-10);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   return ret;
}

bool testAdd3D3()
{
   // Tests the first add method to do scalation of 3D Histograms

   Double_t c1 = r.Rndm();

   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",
                       numberOfBins, minRange, maxRange,
                       numberOfBins+1, minRange, maxRange,
                       numberOfBins+2, minRange, maxRange);

   h1->Sumw2();h2->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);
      Int_t binx = h1->GetXaxis()->FindBin(x);
      Int_t biny = h1->GetYaxis()->FindBin(y);
      Int_t binz = h1->GetZaxis()->FindBin(z);
      Double_t area = h1->GetXaxis()->GetBinWidth( binx ) *
                      h1->GetYaxis()->GetBinWidth( biny ) *
                      h1->GetZaxis()->GetBinWidth( binz );
      h2->Fill(x, y, z, c1 / area);
   }

   TH3D* h3 = new TH3D("t1D1-h3", "h3=c1*h1",
                       numberOfBins, minRange, maxRange,
                       numberOfBins+1, minRange, maxRange,
                       numberOfBins+2, 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
   h2->ResetStats();

   bool ret = equals("Add2D3", h2, h3, cmpOptStats, 1E-10);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   return ret;
}

bool testAdd2D1()
{
   // Tests the first Add method for 2D Histograms

   Double_t c1 = r.Rndm();
   Double_t c2 = r.Rndm();

   TH2D* h1 = new TH2D("t2D1-h1", "h1",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);

   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",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, 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);
      h3->Fill(x, y, 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);
      h2->Fill(x, y,  1.0);
      h3->Fill(x, y, c2);
   }

   TH2D* h4 = new TH2D("t2D1-h4", "h4=c1*h1+c2*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   h4->Add(h1, h2, c1, c2);
   bool ret = equals("Add2D1", h3, h4, cmpOptStats , 1E-10);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testAdd2DProfile1()
{
   // Tests the first Add method for 1D Profiles

   Double_t c1 = r.Rndm();
   Double_t c2 = r.Rndm();

   TProfile2D* p1 = new TProfile2D("t2D1-p1", "p1",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile2D* p2 = new TProfile2D("t2D1-p2", "p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile2D* p3 = new TProfile2D("t2D1-p3", "p3=c1*p1+c2*p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 2, 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 + 2, 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()
{
   // Tests the second Add method for 2D Histograms

   Double_t c2 = r.Rndm();

   TH2D* h1 = new TH2D("t2D2-h1", "h1",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);

   TH2D* h2 = new TH2D("t2D2-h2", "h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);

   TH2D* h3 = new TH2D("t2D2-h3", "h3=h1+c2*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, 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);
      h3->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);
      h3->Fill(x, y, c2);
   }

   h1->Add(h2, c2);
   bool ret = equals("Add2D2", h3, h1, cmpOptStats, 1E-10);
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testAdd2DProfile2()
{
   // Tests the second Add method for 2D Profiles

   Double_t c2 = r.Rndm();

   TProfile2D* p1 = new TProfile2D("t2D2-p1", "p1",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile2D* p2 = new TProfile2D("t2D2-p2", "p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile2D* p3 = new TProfile2D("t2D2-p3", "p3=p1+c2*p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 2, 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()
{
   // Tests the first Add method for 3D Histograms

   Double_t c1 = r.Rndm();
   Double_t c2 = r.Rndm();

   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",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 1, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);

   TH3D* h3 = new TH3D("t3D1-h3", "h3=c1*h1+c2*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 1, minRange, maxRange,
                       numberOfBins + 2, 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);
      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(x, y, z,  1.0);
      h3->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);
      h2->Fill(x, y, z,  1.0);
      h3->Fill(x, y, z, c2);
   }

   TH3D* h4 = new TH3D("t3D1-h4", "h4=c1*h1+c2*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 1, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   h4->Add(h1, h2, c1, c2);
   bool ret = equals("Add3D1", h3, h4, cmpOptStats, 1E-10);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testAdd3DProfile1()
{
   // Tests the second Add method for 3D Profiles

   Double_t c1 = r.Rndm();
   Double_t c2 = r.Rndm();

   TProfile3D* p1 = new TProfile3D("t3D1-p1", "p1",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 1, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile3D* p2 = new TProfile3D("t3D1-p2", "p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 1, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile3D* p3 = new TProfile3D("t3D1-p3", "p3=c1*p1+c2*p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 1, minRange, maxRange,
                                   numberOfBins + 2, 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 + 1, minRange, maxRange,
                                   numberOfBins + 2, 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()
{
   // Tests the second Add method for 3D Histograms

   Double_t c2 = r.Rndm();

   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",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 1, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);

   TH3D* h3 = new TH3D("t3D2-h3", "h3=h1+c2*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 1, minRange, maxRange,
                       numberOfBins + 2, 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);
      Double_t z = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(x, y, z,  1.0);
      h3->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);
      h3->Fill(x, y, z, c2);
   }

   h1->Add(h2, c2);
   bool ret = equals("Add3D2", h3, h1, cmpOptStats, 1E-10);
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testAdd3DProfile2()
{
   // Tests the second Add method for 3D Profiles

   Double_t c2 = r.Rndm();

   TProfile3D* p1 = new TProfile3D("t3D2-p1", "p1",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 1, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile3D* p2 = new TProfile3D("t3D2-p2", "p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 1, minRange, maxRange,
                                   numberOfBins + 2, minRange, maxRange);

   TProfile3D* p3 = new TProfile3D("t3D2-p3", "p3=p1+c2*p2",
                                   numberOfBins, minRange, maxRange,
                                   numberOfBins + 1, minRange, maxRange,
                                   numberOfBins + 2, 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;
}

template<typename HIST>
bool testAddHn()
{
   // Tests the Add method for n-dimensional Histograms

   Double_t c = r.Rndm();

   Int_t bsize[] = { TMath::Nint( r.Uniform(1, 5) ),
                          TMath::Nint( r.Uniform(1, 5) ),
                          TMath::Nint( r.Uniform(1, 5) )};
   Double_t xmin[] = {minRange, minRange, minRange};
   Double_t xmax[] = {maxRange, maxRange, maxRange};

   HIST* s1 = new HIST("tS-s1", "s1", 3, bsize, xmin, xmax);
   HIST* s2 = new HIST("tS-s2", "s2", 3, bsize, xmin, xmax);
   HIST* s3 = new HIST("tS-s3", "s3=s1+c*s2", 3, bsize, xmin, xmax);

   s1->Sumw2();s2->Sumw2();s3->Sumw2();

   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
      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);
      s3->Fill(points);
   }

   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
      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 );
      s2->Fill(points);
      s3->Fill(points, c);
   }

   s1->Add(s2, c);
   bool ret = equals(TString::Format("AddHn<%s>", HIST::Class()->GetName()), s3, s1, cmpOptStats , 1E-10);
   delete s2;
   delete s3;
   return ret;
}

bool testMul1()
{
   // Tests the first Multiply method for 1D Histograms

   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);

   h1->Sumw2();h2->Sumw2();h3->Sumw2();

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value, 1.0);
   }

   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h2->Fill(value,  1.0);
      h3->Fill(value,  c1*c2*h1->GetBinContent( h1->GetXaxis()->FindBin(value) ) );
   }

   // h3 has to be filled again so that the erros are properly calculated
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h3->Fill(value,  c1*c2*h2->GetBinContent( h2->GetXaxis()->FindBin(value) ) );
   }

   // No the bin contents has to be reduced, as it was filled twice!
   for ( Int_t bin = 0; bin <= h3->GetNbinsX() + 1; ++bin ) {
      h3->SetBinContent(bin, h3->GetBinContent(bin) / 2 );
   }

   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);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}