// ----------------------------------------------------------------------
//
// HepHBookHist2D.cc - implementation of HBook two-dimensional histogram
//
// HepHBookHist2D is the implementation of the HBook two-dimensional
// histogram. It is derived from HepHist2D (the generic interface) and
// HepHBookHist. It is mainly a C++ wrapper around the CERNlib HBook
// Fortran histogramming libraries.
//
// No HepHBookHist2D objects are ever directly
// instantiated by the programmer. When creating new HepHist* objects,
// the manager is in charge of creating the appropriate HepHBook object.
//
// To use generate HBook histograms, HepHBookFileManager must be used:
//
//   int lunit = 10;
//   HepFileManager* m = new HepHBookFileManager("filename.rz", lunit);
//    :
//    :
//   HepHist2D* h(m->hist2D("Title", 100, 0.0, 23.0, 100, 0.0, 23.0));
//   float x, y, weight;
//    :
//   [gather data]
//    :
//   h->accumulate(x, y, weight);
//
//
// Note that the constructor can take an ID number which is used by
// the HBook implementation of HepHist2D. If omitted, the ID number is
// generated automatically (which is recommended). The ID number is
// available through the id() method.
//
// Bob Jacobsen
//   30-May-1997  Walter Brown   Added cloning functions
//   29-May-1997  Jason Luther   Added comments
//   04-Jun-1997  Jason Luther   Rewrote comments
//   04-Aug-1998  Philippe Canal Added weight, entries, sum, etc..
//   25-Feb-1999  J. Marraffino  Added increment operators for histograms
//
// ----------------------------------------------------------------------


#ifndef HEPTRACE_H
#include "HepTuple/HepTrace.h"
#endif

#ifndef HEPHBOOKHIST2D_H
#include "HepTuple/HepHBookHist2D.h"
#endif

#ifndef HEPRAWHIST_H
#include "HepTuple/HepRawHist.h"
#endif

#ifndef HEPHBOOKWRAPPERS_H
#include "HepTuple/HepHBookWrappers.h"
#include "ZMtools/pretendToUse.h"
#endif

#include "ZMutility/iostream"

#include <string>
USING( std::string )


ZM_BEGIN_NAMESPACE( zmht )	/*  namespace zmht  {  */

// create a 2D histogram with the given title and binning
HepHBookHist2D::HepHBookHist2D(
  HepHBookFileManager * manager
, const string& title
, const int nBinsX, const float lowX, const float highX
, const int nBinsY, const float lowY, const float highY
, const int idReq
) :
  HepHist2D( manager, title,
             nBinsX, lowX, highX,
             nBinsY, lowY, highY,
             idReq )
{
  HEP_DEBUG( "HepHBookHist2D::constructor( \"" << title << "\", " << *id_ << " )" );
}  //  HepHBookHist2D::HepHBookHist2D()


// Used for dummy construction ...
// the object is then not valid ...
HepHBookHist2D::HepHBookHist2D() : HepHist2D() {
  HEP_DEBUG( "HepHBookHist2D::constructor( )" );
}

void HepHBookHist2D::silenceCompiler() {
  pretendToUse(&(__cf__HEXIST));
  pretendToUse(&(__cf__HSUM));
  pretendToUse(&(__cf__HI));
  pretendToUse(&(__cf__HIE));
  pretendToUse(&(__cf__HIF));
  pretendToUse(&(__cf__HIJ));
  pretendToUse(&(__cf__HIJE));
  pretendToUse(&(__cf__HMIN));
  pretendToUse(&(__cf__HMAX));
  pretendToUse(&(__cf__HRNDM1));
  pretendToUse(&(__cf__HSPFUN));
  pretendToUse(&(__cf__HLNXTF));
  pretendToUse(&(__cf__HSTATI));
  pretendToUse(&(__cf__HX));
  pretendToUse(&(__cf__HXE));
  pretendToUse(&(__cf__HXY));
  pretendToUse(&(__cf__ISCWN));
  pretendToUse(&c2fstrv);
  pretendToUse(&f2cstrv);
  pretendToUse(&vkill_trailing);
  pretendToUse(&num_elem);
}

HepHBookHist2D::~HepHBookHist2D()  {
  HEP_DEBUG( "HepHBookHist2D::destructor( \"" << *title_ << "\", " << *id_ << " )" );
}


void HepHBookHist2D::accumulate (
  const float x
, const float y
, const float weight
)  {

  if ( mayUse() )  {
    _accumulate( x, y, weight );

    HBookPushd( dir() );
    HFILL( *id_, x, y, weight );
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem("attempt to accumulate() to unmanaged histogram."));
    return;
  }
}


float HepHBookHist2D::bin(
  const int binNumX
, const int binNumY
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    float result = HIJ ( *id_, binNumX, binNumY );
    HBookPopd();
    return result;
   }
  else  {
    ZMthrow(ZMxHepUnmanagedItem("attempt to bin() from unmanaged histogram."));
    return 0.0;
  }
}

#define BADCALL(name)							\
  if (! mayUse() ) {							\
    ZMthrow(ZMxHepUnmanagedItem(					\
	    string("attempt to " #name "() an unmanaged histogram ")	\
	    +title()+string(".")));					\
    return 0;								\
  }

int   HepHBookHist2D::entries() const { // Number of entries in bins
  BADCALL(entries);
  int result = 0;

  HBookPushd( dir() );
  HNOENT(*id_, result);
  HBookPopd();

  return result; 
}

float HepHBookHist2D::weight() const {   // Sum of all weights
  BADCALL(weight);
  HBookPushd( dir() );
  float result = HSUM(*id_);
  HBookPopd();
  return result;
}

float HepHBookHist2D::weight2() const { // retrieve sum of weight^2
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to weight2() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  float sum = 0;
  float weight;
  int numBinsX = this->nBinsX();
  int numBinsY = this->nBinsY();
  HBookPushd( dir() );
  for (int i=1;i<numBinsX+1;i++) {
    for (int j=1;j<numBinsY+1;j++) {
    weight = HIJ( *id_, i, j );
    sum += weight * weight;
    }
  }
  HBookPopd();
  return sum;
}

float HepHBookHist2D::sumX() const {  // retrieve weighted sum of X
  BADCALL(sumX);
  float sum = 0;
  float lowX = minX();
  int numBinsX = this->nBinsX();
  int numBinsY = this->nBinsY();
  float deltaX = (maxX()-lowX)/numBinsX;
  HBookPushd( dir() );
  // loop over X
  for (int i=1;i<numBinsX+1;i++) {
    float theX = lowX + deltaX * (i-0.5);
    for (int j=1;j<numBinsY+1;j++)
      sum += theX * HIJ(*id_, i, j);
  } // end for i
  HBookPopd();
  return sum;
}

float HepHBookHist2D::sumY() const {  // retrieve weighted sum of Y 
  BADCALL(sumY);
  float sum = 0;
  float lowY = minY();
  int numBinsX = this->nBinsX();
  int numBinsY = this->nBinsY();
  float deltaY = (maxY()-lowY)/numBinsY;
  HBookPushd( dir() );
  // loop over X
  for (int i=1;i<numBinsY+1;i++) {
    float theY = lowY + deltaY * (i-0.5);
    for (int j=1;j<numBinsX+1;j++)
      sum += theY * HIJ(*id_, j, i);
  } // end for i
  HBookPopd();
  return sum;
}

float HepHBookHist2D::sumX2() const {  // retrieve weighted sum of X^2
  BADCALL(sumX2);
  float sum = 0;
  float lowX = minX();
  int numBinsX = this->nBinsX();
  int numBinsY = this->nBinsY();
  float deltaX = (maxX()-lowX)/numBinsX;
  HBookPushd( dir() );
  // loop over X
  for (int i=1;i<numBinsX+1;i++) {
    float theX = lowX + deltaX * (i-0.5);
    for (int j=1;j<numBinsY+1;j++)
      sum += theX * theX * HIJ(*id_, i, j);
  } // end for i
  HBookPopd();
  return sum;
}

float HepHBookHist2D::sumXY() const {  // retrieve weighted sum of X*Y
  BADCALL(sumXY);
  float sum = 0;
  float lowX = minX();
  float lowY = minY();
  int numBinsX = this->nBinsX();
  int numBinsY = this->nBinsY();
  float deltaX = (maxX()-lowX)/numBinsX;
  float deltaY = (maxY()-lowY)/numBinsY;
  HBookPushd( dir() );

  float* theY = new float[numBinsY+2];
  for (int k=0;k<numBinsY+2;k++)
    theY[k] = lowY + deltaY * (k-0.5);

  // loop over X
  for (int i=1;i<numBinsX+1;i++) {
    float theX = lowX + deltaX * (i-0.5);
    // loop over Y
    for (int j=1;j<numBinsY+1;j++) {
      sum += theY[j] * theX * HIJ(*id_, i, j);
    }
  } // end for i
  HBookPopd();
  return sum;

}

float HepHBookHist2D::sumY2() const {  // retrieve weighted sum of Y^2
  BADCALL(sumY2);
  float sum = 0;
  float lowY = minY();
  int numBinsX = this->nBinsX();
  int numBinsY = this->nBinsY();
  float deltaY = (maxY()-lowY)/numBinsY;
  HBookPushd( dir() );
  // loop over Y
  for (int i=1;i<numBinsY+1;i++) {
    float theY = lowY + deltaY * (i-0.5);
    for (int j=1;j<numBinsX+1;j++)
      sum += theY * theY * HIJ(*id_, j, i);
  } // end for i
  HBookPopd();
  return sum;
}

void HepHBookHist2D::reset()  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    char* chopt = " ";
    HRESET ( *id_, chopt );
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem("attempt to reset() unmanaged histogram."));
    return;
  }
}

float HepHBookHist2D::binError(                // retrieve error of given bin
  const int binNumX
, const int binNumY
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    float result = HIJE( *id_, binNumX, binNumY );
    HBookPopd();
    return result;
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access to unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

void HepHBookHist2D::getContents(             // get values for all in-range bins
  float * data
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );

//  Can't use HUNPAK here. It returns the 2D array of contents in
//  Fortran order and we need the transpose. Do it here, by hand!

    int k;
    int numBinsX = nBinsX();
    int numBinsY = nBinsY();
    for( int i = 0; i<numBinsX; ++i ) {
      for( int j = 0; j<numBinsY; ++j ) {
        k = j + numBinsY*i;
        data[k] = HIJ( *id_, i+1, j+1 );
      }
    }
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHist2D::putContents(             // set values for all in-range bins
  float * data
) const  {
  if ( mayUse() )  {
    HBookPushd( dir() );
    int numBinsX = nBinsX();
    int numBinsY = nBinsY();
    float * tcontents = new float[numBinsX*numBinsY];

// The contents array is presumed to be in C order and we
// need it in Fortran order. Do the transpose here by hand.

    int k;
    int l=0;
    for( int i = 0; i<numBinsX; ++i ) {
      for( int j = 0; j<numBinsY; ++j ) {
        k = i + j*numBinsX;
        tcontents[k] = data[l];
        l++;
      }
    }

// Move the bin contents 

    HPAK ( id(), tcontents );
    delete [] tcontents;
    HBookPopd();
  } else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHist2D::getErrors(               // get errors for all bins
  float * data
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
//  HUNPKE( *id_, *data, "HIST", 1 );

//  Can't use HUNPAKE here. It returns the 2D array of contents in
//  Fortran order and we need the transpose. Do it here, by hand!

    int k;
    int numBinsX = nBinsX();
    int numBinsY = nBinsY();
    for( int i = 0; i<numBinsX; ++i ) {
      for( int j = 0; j<numBinsY; ++j ) {
        k = j + numBinsY*i;
        data[k] = HIJE( *id_, i+1, j+1 );
      }
    }
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHist2D::putErrors(             // set errors for all in-range bins
  float * data
) const  {
  if ( mayUse() )  {
    HBookPushd( dir() );
    int numBinsX = nBinsX();
    int numBinsY = nBinsY();
    float * terrors = new float[numBinsX*numBinsY];

// The errors array is presumed to be in C order and we
// need it in Fortran order. Do the transpose here by hand.

    int k;
    int l=0;
    for( int i = 0; i<numBinsX; ++i ) {
      for( int j = 0; j<numBinsY; ++j ) {
        k = i + j*numBinsX;
        terrors[k] = data[l];
        l++;
      }
    }

// Move the bin contents 

    HPAKE ( id(), terrors );
    delete [] terrors;
    HBookPopd();
  } else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHist2D::getStatistics(		// get histogram-wide statistics
      int & numEntries
 ,  float & sumW
 ,  float & sumW2
 , double & sumWX
 , double & sumWX2
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    numEntries = entries();
    sumW = weight();
    sumW2 = weight2();
    sumWX = (double)sumX();
    sumWX2 = (double)sumX2();
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHist2D::getOverflows(		// get over/under-flow counts
   float * numXOverYOver
 , float * numXOverYUnder
 , float * numXUnderYUnder
 , float * numXUnderYOver
 , float * numTop
 , float * numRight
 , float * numBottom
 , float * numLeft
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    *numXOverYOver   = HIJ( *id_, nBinsX()+1, nBinsY()+1 );
    *numXOverYUnder  = HIJ( *id_, nBinsX()+1,          0 );
    *numXUnderYUnder = HIJ( *id_,          0,          0 );
    *numXUnderYOver  = HIJ( *id_,          0, nBinsY()+1 );

    for( int i = 0; i < nBinsX(); ++i ) {
      numTop[i]    = HIJ( *id_, i+1, nBinsY()+1 );
      numBottom[i] = HIJ( *id_, i+1,          0 );
    }
    for( int i = 0; i < nBinsY(); ++i ) {
      numRight[i] = HIJ( *id_, nBinsX()+1, i+1 );
      numLeft[i]  = HIJ( *id_,          0, i+1 );
    }
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHist2D::setXTitle( const std::string& label )
{  // Dummy function since there is no implementation in HBook for this.
}

void HepHBookHist2D::setYTitle( const std::string& label )
{  // Dummy function since there is no implementation in HBook for this.
}

void HepHBookHist2D::setZTitle( const std::string& label )
{  // Dummy function since there is no implementation in HBook for this.
}

HepRawHist HepHBookHist2D::importHist2D( const HepHist2D & h ) {

// Gather up all the data we need from h

  int numBinsX = h.nBinsX();
  int numBinsY = h.nBinsY();
  float minX = h.minX();
  float maxX = h.maxX();
  float minY = h.minY();
  float maxY = h.maxY();

// Contents and Errors

  float * contents  = new float[numBinsX*numBinsY];
  float * tcontents = new float[numBinsX*numBinsY];
  float * errors    = new float[numBinsX*numBinsY];
  float * terrors   = new float[numBinsX*numBinsY];
  h.getContents( contents );
  h.getErrors( errors );

// The contents and errors arrays are presumed to be in C order and we
// need them in Fortran order. Do the transpose here by hand.

  int k;
  int l=0;
  for( int i = 0; i<numBinsX; ++i ) {
    for( int j = 0; j<numBinsY; ++j ) {
      k = i + j*numBinsX;
      tcontents[k] = contents[l];
      terrors[k] = errors[l];
      l++;
    }
  }
  delete [] contents;
  delete [] errors;

// Overflows and Underflows

  float numXOverYOver;
  float numXOverYUnder;
  float numXUnderYUnder;
  float numXUnderYOver;
  float * numTop    = new float[numBinsX];
  float * numRight  = new float[numBinsY];
  float * numBottom = new float[numBinsX];
  float * numLeft   = new float[numBinsY];

  h.getOverflows( &numXOverYOver, &numXOverYUnder, &numXUnderYUnder, &numXUnderYOver,
                  numTop,        numRight,       numBottom,       numLeft );

// Overall Statistics

  int numEntries;
  float sumW;   
  float sumW2; 
  double sumWX;
  double sumWX2;
  h.getStatistics( numEntries, sumW, sumW2, sumWX, sumWX2 );

// Manufacture a temporary histogram with all the characteristics of h, in the
// same HBook directory as h, but in the same implementation as the result.
// Note that this thing is unmanaged so we need to talk with it directly.

// Find an unused id

  int id;
  HBookPushd( dir() );
  for( id = 1000; HEXIST(id) != 0; ++id) { }

// Book a HBook 2D histogram to park stuff

  char* title = "Temporary hist2D for importHist2D";
  HBOOK2( id, title, numBinsX, minX, maxX,
                     numBinsY, minY, maxY, 0.0 );

// Instantiate a HepRawHist object to pass back to say where the
// Raw Histogram lives and what it's called

  HepRawHist temp( id );

// Move the bin contents and bin errors

  HPAK ( id, tcontents );
  HPAKE( id, terrors );
  delete [] tcontents;
  delete [] terrors;

// Then underflows and overflows

  float tempx;
  float tempy;
  if( numXOverYOver > 0 ) {
    tempx = maxX+1.0;
    tempy = maxY+1.0;
    HFILL( id, tempx, tempy, numXOverYOver );
  }
  if( numXOverYUnder > 0 ) {
    tempx = maxX+1.0;
    tempy = minY-1.0;
    HFILL( id, tempx, tempy, numXOverYUnder );
  }
  if( numXUnderYUnder > 0 ) {
    tempx = minX-1.0;
    tempy = minY-1.0;
    HFILL( id, tempx, tempy, numXUnderYUnder );
  }
  if( numXUnderYOver > 0 ) {
    tempx = minX-1.0;
    tempy = maxY+1.0;
    HFILL( id, tempx, tempy, numXUnderYOver );
  }

  float binWidthX = ( maxX - minX )/(float)numBinsX;
  float binWidthY = ( maxY - minY )/(float)numBinsY;

  for(int i = 0; i < numBinsX; ++i ) {
    tempx = minX + ((float)i+0.5)*binWidthX;
    tempy = maxY+1.0;
    HFILL( id, tempx, tempy, numTop[i] );
    tempy=minY-1.0;
    HFILL( id, tempx, tempy, numBottom[i] );
  }

  for(int i = 0; i < numBinsY; ++i ) {
    tempx = maxX+1.0;
    tempy = minY + ((float)i+0.5)*binWidthY;
    HFILL( id, tempx, tempy, numRight[i] );
    tempx = minX-1.0;
    HFILL( id, tempx, tempy, numLeft[i] );
  }
  delete [] numTop;
  delete [] numBottom;
  delete [] numRight;
  delete [] numLeft;

// Now set the overall statistics

// HBook only wants numEntries for 2D plots. Set sumW2 negative to signal that.
  sumW2 = -1000.0;
  PUTSTAT( temp.hid, numEntries, sumW, sumW2, sumWX, sumWX2 );

// Since we're in control here, let's postpone the Pop and save some time

//HBookPopd();
  return temp;
}

void HepHBookHist2D::addHist2D( const HepRawHist & temp ) {

// If we got here, the temporary histogram object referred to by temp is
// known to be valid. On the other hand, it's worth being careful. 

// Since we postponed the Pop in importHist2D, we shouldn't Push here

//HBookPushd( dir() );
  char* operation = "+";
  if( temp.hid != 0 ) HOPERA( *id_, operation, temp.hid, *id_, 1.0, 1.0 );

// Postpone the Pop again

//HBookPopd ();
}

void HepHBookHist2D::subtractHist2D( const HepRawHist & temp ) {

// If we got here, the temporary histogram object referred to by temp is
// known to be valid. On the other hand, it's worth being careful. 

// Since we postponed the Pop in importHist2D, we shouldn't Push here

//HBookPushd( dir() );
  char* operation = "-";
  if( temp.hid != 0 ) HOPERA( *id_, operation, temp.hid, *id_, 1.0, 1.0 );

// Postpone the Pop again

//HBookPopd ();
}

void HepHBookHist2D::multiplyHist2D( const HepRawHist & temp ) {

// If we got here, the temporary histogram object referred to by temp is
// known to be valid. On the other hand, it's worth being careful. 

// Since we postponed the Pop in importHist2D, we shouldn't Push here

//HBookPushd( dir() );
  char* operation = "*";
  if( temp.hid != 0 ) HOPERA( *id_, operation, temp.hid, *id_, 1.0, 1.0 );

// Postpone the Pop again

//HBookPopd ();
}

void HepHBookHist2D::divideHist2D( const HepRawHist & temp ) {

// If we got here, the temporary histogram object referred to by temp is
// known to be valid. On the other hand, it's worth being careful. 

// Since we postponed the Pop in importHist2D, we shouldn't Push here

//HBookPushd( dir() );
  char* operation = "/";
  if( temp.hid != 0 ) HOPERA( *id_, operation, temp.hid, *id_, 1.0, 1.0 );

// Postpone the Pop again

//HBookPopd ();
}

void HepHBookHist2D::removeImportedHist2D( const HepRawHist & temp ) {

// If we got here, the temporary histogram object referred to by temp is
// known to be valid. On the other hand, it's worth being careful. If we
// pass a zero to HDELET, it will happily (and silently) eat everything!

// Since we postponed the Pop in importHist2D, we shouldn't Push here

//HBookPushd( dir() );
  if( temp.hid != 0 ) HDELET( temp.hid );

// Now, we're really done. Finally, do the Pop we've been postponing.

  HBookPopd ();
}

ZM_END_NAMESPACE( zmht )	/*  }  // namespace zmht  */