#ifndef ZMENVIRONMENT_H
#  include "ZMutility/ZMenvironment.h"
#endif
#include "CLHEP/Random/Randomize.h"
#include "HepTuple/HepRootHist1D.h"
#include "HepTuple/HepRootHist2D.h"
#include "HepTuple/HepRootHistProf.h"
#include "HepTuple/HepRootFileManager.h"
#include "HepTuple/HepHBookHist1D.h"
#include "HepTuple/HepHBookHist2D.h"
#include "HepTuple/HepHBookHistProf.h"
#include "HepTuple/HBookFile.h"
#include "ZMutility/iostream"
#include "ZMutility/iomanip"
#include "ZMutility/cmath"

USING( namespace std )

int main()  {

  RanluxEngine  Ranlux(19735091);
  RandChiSquare Rchi(Ranlux);
  HBookFile * aManager = new HBookFile( "HBtoR.rz" );
  HepRootFileManager * bManager = new HepRootFileManager( "HBtoR.root" );
  int nbins = 50;
  int samples = 5000;
  double dev;

//  Instantiate a gang of 1D histograms.

  HepHist1D & hchi30  = aManager->hist1D("HChi30",      nbins, 0.0f,100.0f);
  HepHist1D & hchi50  = bManager->hist1D("HChi50",      nbins, 0.0f,100.0f);
  HepHist1D & hsum    = bManager->hist1D("HSum",        nbins, 0.0f,100.0f);
  HepHist1D & hdiff   = bManager->hist1D("HDifference", nbins, 0.0f,100.0f);
  HepHist1D & hprod   = bManager->hist1D("HProduct",    nbins, 0.0f,100.0f);
  HepHist1D & hratio  = bManager->hist1D("HRatio",      nbins, 0.0f,100.0f);
  HepHist1D & hcopy   = bManager->hist1D("HCopy",       nbins, 0.0f,100.0f);
  HepHist1D & hflip   = bManager->hist1D("HFlip",       nbins, 0.0f,100.0f);

//  Make some entries.  Fill with random numbers from various distributions.

  cout << "\nTests for 1D histograms" << endl;
  for( int j=0; j<samples; j++ ) {
    dev = Rchi.fire( 30.0 );
    hchi30.accumulate( (float)dev );
  }
  for( int j=0; j<samples; j++ ) {
    dev = Rchi.fire( 50.0 );
    hchi50.accumulate( (float)dev );
    dev = Rchi.fire( 50.0 );
    hchi50.accumulate( (float)dev );
  }

// Throw in some overflows and underflows

  for( int i = 0; i<3; i++ ) {
    hchi30.accumulate( 150.0f );
    hchi50.accumulate( -2.0f );
    hchi50.accumulate( -4.0f );
  }

// Now form the sum of chi30 and chi50

  cout << "Testing +=" << endl;
  hsum += hchi30;
  hsum += hchi50;

// And the difference of chi30 and chi50

  cout << "Testing -=" << endl;
  hdiff += hchi30;
  hdiff -= hchi50;

// And the product of chi30 and chi50

  cout << "Testing *=" << endl;
  hprod += hchi30;
  hprod *= hchi50;

// And the ratio of chi30 and chi50

  cout << "Testing /=" << endl;
  hratio += hchi30;
  hratio /= hchi50;

// and finally some test things

  cout << "Poor man's copy" << endl;
  hcopy += hchi30;
  cout << "Poor man's copy with negate" << endl;
  hflip -= hchi50;

// ----------------------------------------------------------------------

  int nbinsx = 20;
  int nbinsy = 10;
  float x1[16] = {  5.5,  5.5, -5.5, -5.5, 25.0, 25.0,-25.0,-25.0,
                   -7.5,  7.5, -3.5,  3.5, 25.0, 25.0,-25.0,-25.0 };
  float y1[16] = {  5.0, -5.0,  5.0, -5.0, -6.0,  6.0, -4.0,  4.0,
                   25.0, 25.0,-25.0,-25.0, 25.0,-25.0,-25.0, 25.0 };
  int mult1[16]= {    2,    4,    6,    8,    1,    2,    3,    4,
                      5,    6,    7,    8,    1,    2,    3,    4 };
  float x2[16] = {  5.5,  5.5, -5.5, -5.5,-25.0,-25.0, 25.0, 25.0,
                    7.5, -7.5,  3.5, -3.5,-25.0,-25.0, 25.0, 25.0 };
  float y2[16] = {  5.0, -5.0,  5.0, -5.0,  6.0, -6.0,  4.0, -4.0,
                  -25.0,-25.0, 25.0, 25.0,-25.0, 25.0, 25.0,-25.0 };
  int mult2[16]= {    1,    2,    3,    4,    4,    3,    2,    1,
                      8,    7,    6,    5,    4,    3,    2,    1 };

//  Instantiate a gang of 2D histograms.

  HepHist2D & lego1 = aManager->hist2D("Slego1",     nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & lego2 = bManager->hist2D("Slego2",     nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & ssum  = bManager->hist2D("SSum",       nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & sdiff = bManager->hist2D("SDifference",nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & sprod = bManager->hist2D("SProduct",   nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & sratio= bManager->hist2D("SRatio",     nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & scopy = bManager->hist2D("SCopy",      nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);
  HepHist2D & sflip = bManager->hist2D("SFlip",      nbinsx, -10.0f,10.0f,
                                                     nbinsy, -10.0f,10.0f,0);

//  Make some entries.  Fill with random numbers from various distributions.

  cout << "\nTests for 2D histograms" << endl;
  float w = 1.1;
  for( int i=0; i<16; i++ ) {
    for( int j=0; j<mult1[i]; j++ ) {
      lego1.accumulate( x1[i], y1[i], w );
    }
    for( int j=0; j<mult2[i]; j++ ) {
      lego2.accumulate( x2[i], y2[i], w );
    }
  }

// Now form the sum of lego1 and lego2

  cout << "Testing +=" << endl;
  ssum += lego1;
  ssum += lego2;

// And the difference of lego1 and lego2

  cout << "Testing -=" << endl;
  sdiff += lego1;
  sdiff -= lego2;

// And the product of lego1 and lego2

  cout << "Testing *=" << endl;
  sprod += lego1;
  sprod *= lego2;

// And the ratio of lego1 and lego2

  cout << "Testing /=" << endl;
  sratio += lego1;
  sratio /= lego2;

// and finally some test things

  cout << "Poor man's copy" << endl;
  scopy += lego1;
  cout << "Poor man's copy with negate" << endl;
  sflip -= lego2;

// ----------------------------------------------------------------------

//  Instantiate a gang of Profile histograms.

  HepHistProf & pchi1 = aManager->histProf("PChi1",       nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & pchi2 = bManager->histProf("PChi2",       nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & psum  = bManager->histProf("PSum",        nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & pdiff = bManager->histProf("PDifference", nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & pprod = bManager->histProf("PProduct",    nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & pratio= bManager->histProf("PRatio",      nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & pcopy = bManager->histProf("PCopy",       nbins, 0.0f,50.0f, -100.0f, 100.0f );
  HepHistProf & pflip = bManager->histProf("PFlip",       nbins, 0.0f,50.0f, -100.0f, 100.0f );

//  Make some entries.  Fill with random numbers from various distributions.

  cout << "\nTests for profile histograms" << endl;
  double x, y;
  samples = 500;
  for( int i=0; i<nbins; ++i ) {
    x = (double)i;
    for( int j=0; j<samples; j++ ) {
      y = Rchi.fire( x+1.0 );
      pchi1.accumulate( (float)x, (float)y );
    }
  }
  for( int i=0; i<nbins-10; ++i ) {
    x = (double)i;
    for( int j=0; j<samples; j++ ) {
      y = Rchi.fire( 50.0-x );
      pchi2.accumulate( (float)x, (float)y );
      y = Rchi.fire(50.0-x );
      pchi2.accumulate( (float)x, (float)y );
    }
  }

// Now form the sum of chi1 and chi2

  cout << "Testing +=" << endl;
  psum += pchi1;
  psum += pchi2;

// And the difference of chi1 and chi2

  cout << "Testing -=" << endl;
  pdiff += pchi1;
  pdiff -= pchi2;

// And the product of chi1 and chi2

  cout << "Testing *=" << endl;
  cout << "Warning. Root has only a dummy Multiply method for profile histograms." << endl;
  pprod += pchi1;
//pprod *= pchi2;

// And the ratio of chi1 and chi2

  cout << "Testing /=" << endl;
  pratio += pchi1;
  pratio /= pchi2;

// and finally some test things

  cout << "Poor man's copy" << endl;
  pcopy += pchi1;
  cout << "Poor man's copy with negate" << endl;
  pflip -= pchi2;

// ----------------------------------------------------------------------

//  Write out the histograms and bail out

  cout << "\nHistograms filled.  Write the file." << endl;
  cout << "Now go look at HBtoR.root with a Root browser." << endl;
  delete aManager;
  delete bManager;

  return 0;
}