stressHistogram.cxx

   return ret;
}

bool testMulVar1()
{
   // Tests the first Multiply 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("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();

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

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

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

   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,  h1->GetBinContent( h1->GetXaxis()->FindBin(value) ) );
   }

   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,  h2->GetBinContent( h2->GetXaxis()->FindBin(value) ) );
   }

   for ( Int_t bin = 0; bin <= h3->GetNbinsX() + 1; ++bin ) {
      h3->SetBinContent(bin, h3->GetBinContent(bin) / 2 );
   }

   h1->Multiply(h2);

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

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

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

   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,  h1->GetBinContent( h1->GetXaxis()->FindBin(value) ) );
   }

   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,  h2->GetBinContent( h2->GetXaxis()->FindBin(value) ) );
   }

   for ( Int_t bin = 0; bin <= h3->GetNbinsX() + 1; ++bin ) {
      h3->SetBinContent(bin, h3->GetBinContent(bin) / 2 );
   }

   h1->Multiply(h2);

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

bool testMul2D1()
{
   // Tests the first Multiply method for 2D Histograms

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

   TH2D* h1 = new TH2D("m2D1_h1", "h1-Title",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   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",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, 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 * 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);
   }

   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,  c1*c2*h1->GetBinContent( h1->GetXaxis()->FindBin(x),
                                               h1->GetYaxis()->FindBin(y) ) );
   }

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

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

   TH2D* h4 = new TH2D("m2D1_h4", "h4=h1*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   h4->Multiply(h1, h2, c1, c2);

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

bool testMul2D2()
{
   // Tests the second Multiply method for 2D Histograms

   TH2D* h1 = new TH2D("m2D2_h1", "h1-Title",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   TH2D* h2 = new TH2D("m2D2_h2", "h2-Title",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   TH2D* h3 = new TH2D("m2D2_h3", "h3=h1*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 2, 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 * 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);
   }

   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,  h1->GetBinContent( h1->GetXaxis()->FindBin(x),
                                         h1->GetYaxis()->FindBin(y) ) );
   }

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

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

   h1->Multiply(h2);

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

bool testMul3D1()
{
   // Tests the first Multiply method for 3D Histograms

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

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

   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,  c1*c2*h1->GetBinContent( h1->GetXaxis()->FindBin(x),
                                                  h1->GetYaxis()->FindBin(y),
                                                  h1->GetZaxis()->FindBin(z) ) );
   }

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

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

   TH3D* h4 = new TH3D("m3D1_h4", "h4=h1*h2",
                       numberOfBins, minRange, maxRange,
                       numberOfBins + 1, minRange, maxRange,
                       numberOfBins + 2, minRange, maxRange);
   h4->Multiply(h1, h2, c1, c2);

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

bool testMul3D2()
{
   // Tests the second Multiply method for 3D Histograms

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

   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, h1->GetBinContent( h1->GetXaxis()->FindBin(x),
                                           h1->GetYaxis()->FindBin(y),
                                           h1->GetZaxis()->FindBin(z) ) );
   }

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

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

   h1->Multiply(h2);

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

template<typename HIST>
bool testMulHn()
{
  // Tests the Multiply method for Sparse Histograms

   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("m3D2-s1", "s1-Title", 3, bsize, xmin, xmax);
   HIST* s2 = new HIST("m3D2-s2", "s2-Title", 3, bsize, xmin, xmax);
   HIST* s3 = new HIST("m3D2-s3", "s3=s1*s2", 3, bsize, xmin, xmax);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   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, 1.0);
   }

   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, 1.0);
      Int_t points_s1[3];
      points_s1[0] = s1->GetAxis(0)->FindBin( points[0] );
      points_s1[1] = s1->GetAxis(1)->FindBin( points[1] );
      points_s1[2] = s1->GetAxis(2)->FindBin( points[2] );
      s3->Fill(points, s1->GetBinContent( points_s1 ) );
   }

   // s3 has to be filled again so that the errors are properly calculated
   r.SetSeed(seed);
   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 );
      Int_t points_s2[3];
      points_s2[0] = s2->GetAxis(0)->FindBin( points[0] );
      points_s2[1] = s2->GetAxis(1)->FindBin( points[1] );
      points_s2[2] = s2->GetAxis(2)->FindBin( points[2] );
      s3->Fill(points, s2->GetBinContent( points_s2 ) );
   }

   // No the bin contents has to be reduced, as it was filled twice!
   for ( Long64_t i = 0; i < s3->GetNbins(); ++i ) {
      Int_t bin[3];
      Double_t v = s3->GetBinContent(i, bin);
      s3->SetBinContent( bin, v / 2 );
   }

   s1->Multiply(s2);

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

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

   TF1* f = new TF1("sin", "sin(x)", minRange - 2, maxRange + 2);

   h1->Sumw2();h2->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);
      h2->Fill(value, f->Eval( h2->GetBinCenter( h2->FindBin(value) ) ) * c1 );
   }

   h1->Multiply(f, c1);

   // stats fails because of the error precision
   int status = equals("MULF H1D", h1, h2); //,cmpOptStats | cmpOptDebug);
   if (cleanHistos) delete h1;
   delete f;
   return status;
}

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

   TF2* f = new TF2("sin2", "sin(x)*cos(y)",
                    minRange - 2, maxRange + 2,
                    minRange - 2, maxRange + 2);
   h1->Sumw2();h2->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);
      h2->Fill(value, f->Eval( h2->GetXaxis()->GetBinCenter( h2->GetXaxis()->FindBin(value) ),
                               h2->GetYaxis()->GetBinCenter( h2->GetYaxis()->FindBin(double(0)) ) )
               * c1 );
   }

   h1->Multiply(f, c1);

   // stats fails because of the error precision
   int status = equals("MULF H1D2", h1, h2); //,cmpOptStats | cmpOptDebug);
   if (cleanHistos) delete h1;
   delete f;
   return status;
}

bool testMulF2D()
{
   Double_t c1 = r.Rndm();

   TH2D* h1 = new TH2D("mf2D_h1", "h1-Title",
                       numberOfBins, minRange, maxRange,
                       numberOfBins, minRange, maxRange);
   TH2D* h2 = new TH2D("mf2D_h2", "h2=h1*c1*f1",
                       numberOfBins, minRange, maxRange,
                       numberOfBins, minRange, maxRange);

   TF1* f = new TF1("sin", "sin(x)", minRange - 2, maxRange + 2);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   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);
      h2->Fill(x, y, f->Eval( h2->GetXaxis()->GetBinCenter( h2->GetXaxis()->FindBin(x) ) ) * c1 );
   }

   h1->Multiply(f, c1);

   // stats fails because of the error precision
   int status = equals("MULF H2D", h1, h2); //, cmpOptStats | cmpOptDebug);
   if (cleanHistos) delete h1;
   delete f;
   return status;
}

bool testMulF2D2()
{
   Double_t c1 = r.Rndm();

   TH2D* h1 = new TH2D("mf2D2_h1", "h1-Title",
                       numberOfBins, minRange, maxRange,
                       numberOfBins, minRange, maxRange);
   TH2D* h2 = new TH2D("mf2D2_h2", "h2=h1*c1*f1",
                       numberOfBins, minRange, maxRange,
                       numberOfBins, minRange, maxRange);

   TF2* f = new TF2("sin2", "sin(x)*cos(y)",
                    minRange - 2, maxRange + 2,
                    minRange - 2, maxRange + 2);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   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);
      h2->Fill(x, y, f->Eval( h2->GetXaxis()->GetBinCenter( h2->GetXaxis()->FindBin(x) ),
                              h2->GetYaxis()->GetBinCenter( h2->GetYaxis()->FindBin(y) ) )
               * c1 );
   }

   h1->Multiply(f, c1);

   // stats fails because of the error precision
   int status = equals("MULF H2D2", h1, h2); //, cmpOptStats | cmpOptDebug);
   if (cleanHistos) delete h1;
   delete f;
   return status;
}

bool testMulF3D()
{
   Double_t c1 = r.Rndm();

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

   TF1* f = new TF1("sin", "sin(x)", minRange - 2, maxRange + 2);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   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);
      h2->Fill(x, y, z, f->Eval( h2->GetXaxis()->GetBinCenter( h2->GetXaxis()->FindBin(x) ) ) * c1 );
   }

   h1->Multiply(f, c1);

   // stats fails because of the error precision
   int status = equals("MULF H3D", h1, h2); //, cmpOptStats | cmpOptDebug);
   if (cleanHistos) delete h1;
   delete f;
   return status;
}

bool testMulF3D2()
{
   Double_t c1 = r.Rndm();

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

   TF2* f = new TF2("sin2", "sin(x)*cos(y)",
                    minRange - 2, maxRange + 2,
                    minRange - 2, maxRange + 2);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   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);
      h2->Fill(x, y, z, f->Eval( h2->GetXaxis()->GetBinCenter( h2->GetXaxis()->FindBin(x) ),
                                 h2->GetYaxis()->GetBinCenter( h2->GetYaxis()->FindBin(y) ) )
               * c1 );
   }

   h1->Multiply(f, c1);

   // stats fails because of the error precision
   int status = equals("MULF H3D2", h1, h2); //, cmpOptStats | cmpOptDebug);
   if (cleanHistos) delete h1;
   delete f;
   return status;
}

template <typename HIST>
bool testMulFND()
{
   const UInt_t nDims = 3;
   Double_t c1 = 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("mfND-s1", "s1-Title", nDims, bsize, xmin, xmax);
   HIST* s2 = new HIST("mfND-s2", "s2=f*s2",  nDims, bsize, xmin, xmax);

   TF1* f = new TF1("sin", "sin(x)", minRange - 2, maxRange + 2);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
      Double_t points[nDims];
      for ( UInt_t i = 0; i < nDims; ++ i )
         points[i] = r.Uniform( minRange * .9 , maxRange * 1.1 );
      s1->Fill(points, 1.0);
      s2->Fill(points, f->Eval( s2->GetAxis(0)->GetBinCenter( s2->GetAxis(0)->FindBin(points[0]) ) ) * c1);
   }

   s1->Multiply(f, c1);

   int status = equals(TString::Format("MULF HND<%s>", HIST::Class()->GetName()), s1, s2);
   delete s1;
   delete f;
   return status;
}

template<typename HIST>
bool testMulFND2()
{
   const UInt_t nDims = 3;
   Double_t c1 = 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("mfND-s1", "s1-Title", nDims, bsize, xmin, xmax);
   HIST* s2 = new HIST("mfND-s2", "s2=f*s2",  nDims, bsize, xmin, xmax);

   TF2* f = new TF2("sin2", "sin(x)*cos(y)",
                    minRange - 2, maxRange + 2,
                    minRange - 2, maxRange + 2);

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents * nEvents; ++e ) {
      Double_t points[nDims];
      for ( UInt_t i = 0; i < nDims; ++ i )
         points[i] = r.Uniform( minRange * .9 , maxRange * 1.1 );
      s1->Fill(points, 1.0);
      s2->Fill(points, f->Eval( s2->GetAxis(0)->GetBinCenter( s2->GetAxis(0)->FindBin(points[0]) ),
                                s2->GetAxis(1)->GetBinCenter( s2->GetAxis(1)->FindBin(points[1]) ) )
                                * c1);
   }

   s1->Multiply(f, c1);

   int status = equals(TString::Format("MULF HND2<%s>", HIST::Class()->GetName()), s1, s2);
   delete s1;
   delete f;
   return status;
}

bool testDivide1()
{
   // Tests the first Divide method for 1D Histograms

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

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value;
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value, 1.0);
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h2->Fill(value,  1.0);
   }
   // avoid bins in h2 with zero content
   for (int i = 0; i < h2->GetSize(); ++i)
      if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);


   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);
   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);
      error -= (2*(c2*c2)/(c1*c1)) * h3->GetBinContent(bin)*h3->GetBinContent(bin)*h2->GetBinError(bin)*h2->GetBinError(bin);
      h4->SetBinError( bin, sqrt(error) );
   }
   h4->ResetStats();
   h1->ResetStats();

   bool ret = equals("Divide1D1", h1, h4, cmpOptStats );
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testDivideVar1()
{
   // Tests the first Divide method for 1D Histograms with variable bin size

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

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

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value;
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value, 1.0);
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h2->Fill(value,  1.0);
   }
   // avoid bins in h2 with zero content
   for (int i = 0; i < h2->GetSize(); ++i)
      if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);


   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);
   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);
      error -= (2*(c2*c2)/(c1*c1)) * h3->GetBinContent(bin)*h3->GetBinContent(bin)*h2->GetBinError(bin)*h2->GetBinError(bin);
      h4->SetBinError( bin, sqrt(error) );
   }
   h4->ResetStats();
   h1->ResetStats();

   bool ret = equals("DivideVar1D1", h1, h4, cmpOptStats);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}


bool testDivideProf1()
{
   // Tests the first Divide method for 1D Profiles

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

   TProfile* p1 = new TProfile("d1D1_p1", "p1-Title", numberOfBins, minRange, maxRange);
   TProfile* p2 = new TProfile("d1D1_p2", "p2-Title", numberOfBins, minRange, maxRange);

   p1->Sumw2();p2->Sumw2();

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


   TProfile* p3 = new TProfile("d1D1_p3", "p3=(c1*p1)/(c2*p2)", numberOfBins, minRange, maxRange);
   p3->Divide(p1, p2, c1, c2);

   // There is no Multiply method to tests. And the errors are wrongly
   // calculated in the TProfile::Division method, so there is no
   // point to make the tests. Once the method is fixed, the tests
   // will be finished.

   return 0;
}

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

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value;
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value, 1.0);
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h2->Fill(value,  1.0);
   }
   // avoid bins in h2 with zero content
   for (int i = 0; i < h2->GetSize(); ++i)
      if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);

   TH1D* h3 = static_cast<TH1D*>( h1->Clone() );
   h3->Divide(h2);

   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);
      error -= 2 * h3->GetBinContent(bin)*h3->GetBinContent(bin)*h2->GetBinError(bin)*h2->GetBinError(bin);
      h4->SetBinError( bin, sqrt(error) );
   }

   h4->ResetStats();
   h1->ResetStats();

   bool ret = equals("Divide1D2", h1, h4, cmpOptStats);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

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

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

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

   UInt_t seed = r.GetSeed();
   // For possible problems
   r.SetSeed(seed);
   for ( Int_t e = 0; e < nEvents; ++e ) {
      Double_t value;
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h1->Fill(value, 1.0);
      value = r.Uniform(0.9 * minRange, 1.1 * maxRange);
      h2->Fill(value,  1.0);
   }
   // avoid bins in h2 with zero content
   for (int i = 0; i < h2->GetSize(); ++i)
      if (h2->GetBinContent(i) == 0) h2->SetBinContent(i,1);

   TH1D* h3 = static_cast<TH1D*>( h1->Clone() );
   h3->Divide(h2);

   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);
      error -= 2 * h3->GetBinContent(bin)*h3->GetBinContent(bin)*h2->GetBinError(bin)*h2->GetBinError(bin);
      h4->SetBinError( bin, sqrt(error) );
   }

   h4->ResetStats();
   h1->ResetStats();

   bool ret = equals("DivideVar1D2", h1, h4, cmpOptStats);
   if (cleanHistos) delete h1;
   if (cleanHistos) delete h2;
   if (cleanHistos) delete h3;
   return ret;
}

bool testDivide2D1()
{
   // Tests the first Divide method for 2D Histograms

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

   TH2D* h1 = new TH2D("d2D1_h1", "h1-Title",