// ----------------------------------------------------------------------
//
// HepHistoHist2D.cc - implementation of Histo two-dimensional histogram
//
// HepHistoHist2D is the implementation of the Histo two-dimensional
// histogram. It is derived from HepHist2D (the generic interface) and
// HepHistoHist. It is mainly a C++ wrapper around the CERNlib Histo
// Fortran histogramming libraries.
//
// No HepHistoHist2D objects are ever directly
// instantiated by the programmer. When creating new HepHist* objects,
// the manager is in charge of creating the appropriate HepHisto object.
//
// To use generate Histo histograms, HepHistoFileManager must be used:
//
//   int lunit = 10;
//   HepFileManager* m = new HepHistoFileManager("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 Histo implementation of HepHist2D. If omitted, the ID number is
// generated automatically (which is recommended). The ID number is
// available through the id() method.
//
// The code presumes that Fortran is callable from C by the (mostly)
// standard BSD convention: all lowercase symbol names with an
// underscore ("_") appended.
//
// 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..
//   21-Feb-2000  Walter Brown   Improved C++ standard compliance
//
// ----------------------------------------------------------------------


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

#ifndef HEPHISTOHIST2D_H
#include "HepTuple/HepHistoHist2D.h"
#endif

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

#include "histoscope.h"

#include "ZMutility/iostream"

#include "ZMutility/cmath"
USING( std::abs )
USING( std::sqrt )

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

// for sprintf
#include <cstdio>

#define HS_MAX_CATEGORY_LENGTH 255
#define HS_MAX_TITLE_LENGTH 80
#define HS_MAX_NAME_LENGTH 80

#include "HepTuple/HepHistoNew.h"

ZM_BEGIN_NAMESPACE( zmht )	/*  namespace zmht  {  */

// create a 2D histogram with the given title and binning
HepHistoHist2D::HepHistoHist2D(
  HepHistoFileManager * 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 ),
  values_(new float[nBinsX*nBinsY])
{
  HEP_DEBUG( "HepHistoHist2D::constructor( \"" << title << "\", " << id_ << " )" );
  string name = title;
  histo_id_ = hs_create_2d_hist( *id_, (char*) title.c_str(),
			(char*) HepHistoFileManager::category(*dir_).c_str(),
				(char*) (name+" x").c_str(),
				(char*) (name+" y").c_str(),
				(char*) (name+" z").c_str(),
				nBinsX, nBinsY,
				lowX, highX,
				lowY, highY);
  if ( histo_id_ == 0 ) {
    ZMthrow(ZMxHepCantMakeHist(string("Histo failed to create the hist 2d: ") 
			       + title ) );
    isValid(false);
  }
}  //  HepHistoHist2D::HepHistoHist2D()

// create a 2D histogram with the given title and binning
HepHistoHist2D::HepHistoHist2D(
  HepHistoFileManager * 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
, const int histo_id
) :
  HepHist2D( manager, title,
             nBinsX, lowX, highX,
             nBinsY, lowY, highY,
             idReq ),
  histo_id_(histo_id),
  values_(new float[nBinsX*nBinsY])
{
  HEP_DEBUG( "HepHistoHist2D::constructor( \"" << title << "\", " << id_ << " )" );
}  //  HepHistoHist2D::HepHistoHist2D()


// clone an existing 2D histogram
HepHistoHist2D::HepHistoHist2D(
  HepHist2D * original
, const string& title
, const int     idReq
) : HepHist2D( original->manager(), title,
               original->nBinsX(), original->minX(), original->maxX(),
               original->nBinsY(), original->minY(), original->maxY(),
               idReq ),
    values_( new float[original->nBinsX()*original->nBinsY()] )
{

  // make the Histoscope item
  int oldID = ((HepHistoHist2D*)original)->histo_id_;
  histo_id_ = hs_copy(oldID, idReq, (char*)title.c_str(), 
		 (char*)HepHistoFileManager::category(original->dir()).c_str() );

}  //  HepHistoHist2D::HepHistoHist2D()


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

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


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

  if ( mayUse() )  {
    _accumulate( x, y, weight );
    hs_fill_2d_hist( histo_id_, x, y, weight );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem("attempt to accumulate() to unmanaged histogram."));
    return;
  }
}

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

  if ( mayUse() )  {
    if ( binNumX <= 0 || binNumX > nBinsX() ||
	 binNumY <= 0 || binNumY > nBinsY() ) {
      float overflows[3][3];
      hs_2d_hist_overflows( histo_id_, overflows );
      int x = 1;
      int y = 1;
      int factor = 1;
      
      if (binNumX <= 0 ) {
    	x = 0;
      } else if (binNumX > nBinsX() ) {
    	x = 2;
      } else {
	factor = nBinsX();
      }
      if (binNumY <= 0) {
    	y = 2;		/* yes, really 2 */
      } else if (binNumY > nBinsY() ) {
    	y = 0;		/* yes, really 0 */
      } else {
	factor = nBinsY();
      }
      return overflows[x][y] / factor;
    } else {
      return hs_2d_hist_bin_value( histo_id_, binNumX-1, binNumY-1 );
    }
  }
  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;								\
  }

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

  if ( mayUse() )  {
    //    int numBins = nBins();
    float result = 0.0;
    /*float * errors = new float[numBins];
    int status = hs_2d_hist_errors( histo_id_, errors, 0 );
    switch(status) {
    case HS_NO_ERRORS:
      result = sqrt(abs(bin(binNum)));	// return best guess
      break;
    case HS_POS_ERRORS:
      result = errors[binNum-1];    
      break;
    case HS_BOTH_ERRORS:
      result = errors[binNum-1];    
      break;
    default:
      char temp[16];
      sprintf(temp,"%d",type);
      ZMthrow(ZMxHepUnknownHist(
			    string("binError: unknown histogram (")
			    +string(temp)+string(")")
			    +string(" in \"")+manager()->fileName()+string("\")")) );
      result = 0.0;
    } // end of switch
    delete [] errors;
    */
    return result;
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

void HepHistoHist2D::getContents(             // get values for all bins
  float * data
) const  {

  if ( mayUse() )  {
    hs_2d_hist_bin_contents( histo_id_, data );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

  if ( mayUse() )  {
    /*  int numBins = nBins();
    float * errors = new float[numBins];
    int status = hs_2d_hist_errors( histo_id_, errors, 0 );
    switch(status) {
    case HS_NO_ERRORS:
      for( int i=0; i<numBins; ++i) {
        data[i] = sqrt(abs(bin(i)));	// return best guess
      }
      break;
    case HS_POS_ERRORS:
      for( int i=0; i<numBins; ++i) {
        data[i] = errors[i];    
      }
      break;
    case HS_BOTH_ERRORS:
      for( int i=0; i<numBins; ++i) {
        data[i] = errors[i];    
      }
      break;
    default:
      char temp[16];
      sprintf(temp,"%d",type);
      ZMthrow(ZMxHepUnknownHist(
			    string("binError: unknown histogram (")
			    +string(temp)+string(")")
			    +string(" in \"")+manager()->fileName()+string("\")")));
      for( int i=0; i<numBins; ++i) {
        data[i] = 0.0;
      }
    } // end of switch
    delete [] errors;
    */
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

/*  New stuff goes here
 *      getStatistics
 *      getOverflows
 *	importHist2D
 *	addHist2D
 *	subtractHist2D
 *	multiplyHist2D
 *	divideHist2D
 *	removeImportedHist2D
 */

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

  numEntries = 0;
  sumW   = 0.0;
  sumW2  = 0.0;
  sumWX  = 0.0;
  sumWX2 = 0.0;

  if ( mayUse() )  {
    numEntries = hs_num_entries( histo_id_ );
    int numBinsX = nBinsX();
    int numBinsY = nBinsY();
    float * contents = new float[numBinsX*numBinsY];
    getContents( contents );

    float thisBin;
    float thisX;
    float xMin = minX();
    float xMax = maxX();
    float binWidthX = ( xMax - xMin )/(float)numBinsX;
    int m = 0;
    for( int i = 0; i < numBinsX; ++i ) {
      for( int k = 0; k < numBinsY; ++k ) {
        m = i*numBinsY + k;
        thisBin = contents[m];
        thisX = xMin + binWidthX*((float)i+0.5);
        sumW   += thisBin;
        sumW2  += thisBin*thisBin;
        sumWX  += thisBin*thisX;
        sumWX2 += thisBin*thisX*thisX;
      }
    }
    delete [] contents;
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

  if ( mayUse() )  {
    float hsover[3][3];
    hs_2d_hist_overflows( histo_id_, hsover );
    *numXUnderYOver  = hsover[0][0];
    *numXOverYOver   = hsover[2][0];
    *numXUnderYUnder = hsover[0][2];
    *numXOverYUnder  = hsover[2][2];

    int numBinsX = this->nBinsX();
    int numBinsY = this->nBinsY();
    float numTopBar    = hsover[1][0]/(float)numBinsX;
    float numBottomBar = hsover[1][2]/(float)numBinsX;
    float numLeftBar   = hsover[0][1]/(float)numBinsY;
    float numRightBar  = hsover[2][1]/(float)numBinsY;
    for( int i = 0; i < numBinsX; ++i ) {
      numTop[i]    = numTopBar;;
      numBottom[i] = numBottomBar;
    }
    for( int i = 0; i < numBinsY; ++i ) {
      numLeft[i]  = numLeftBar;
      numRight[i] = numRightBar;
    }
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

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

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

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

// Gather up all the data we need from h

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

  float * contents  = new float[numBinsX*numBinsY];
  float * errors    = new float[numBinsX*numBinsY];

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

// Contents and errors

  h.getContents( contents );
  h.getErrors( errors );

// 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, but
// in the same implementation as the result.
// Note that this thing is unmanaged so we need to talk with it directly.

// Create a Histoscope 2D histogram to park stuff. Find an unused category
// name to avoid clashes.

  char category[16];
  int * ids = new int[10];
  unsigned int numCategory = 1000; 
  int numMatches = 1;
  while( numMatches != 0 ) {
    sprintf(category,"%X",numCategory++);
    numMatches = hs_list_items( "", category, ids, 10, 1 );
  }
  delete [] ids;

  int uid = 0;
  int tmp_id = 0;
  tmp_id = hs_create_2d_hist( uid, "Temporary hist2D for importHist2D",
             category, "X", "Y", "Z", numBinsX, numBinsY,
              h.minX(), h.maxX(), h.minY(), h.maxY() );

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

  HepRawHist temp( tmp_id );

// Set up to move histogram data.

// Move the bin contents and bin errors

  hs_2d_hist_block_fill( tmp_id, contents, errors, 0 );
  delete [] contents;
  delete [] errors;

// Then underflows and overflows

  if( numXOverYOver > 0 )   hs_fill_2d_hist( tmp_id, xMax+1.0, yMax+1.0,
                                                        numXOverYOver );
  if( numXOverYUnder > 0 )  hs_fill_2d_hist( tmp_id, xMax+1.0, yMin-1.0,
                                                       numXOverYUnder );
  if( numXUnderYUnder > 0 ) hs_fill_2d_hist( tmp_id, xMin-1.0, yMin-1.0,
                                                      numXUnderYUnder );
  if( numXUnderYOver > 0 )  hs_fill_2d_hist( tmp_id, xMin-1.0, yMax+1.0,
                                                       numXUnderYOver );

  float binWidthX = ( xMax - xMin )/(float)numBinsX;
  float binWidthY = ( yMax - yMin )/(float)numBinsY;
  float thisX, thisY;

  for(int i = 0; i < numBinsX; ++i ) {
    thisX = xMin + ((float)i+0.5)*binWidthX;
    thisY = yMax+1.0;
    if( numTop[i] > 0.0 ) hs_fill_2d_hist( tmp_id, thisX, thisY, numTop[i] );
    thisY = yMin-1.0;
    if( numBottom[i] > 0.0 ) hs_fill_2d_hist( tmp_id, thisX, thisY, numBottom[i] );
  }

  for(int i = 0; i < numBinsY; ++i ) {
    thisX = xMax+1.0;
    thisY = yMin + ((float)i+0.5)*binWidthY;
    if( numRight[i] > 0.0 ) hs_fill_2d_hist( tmp_id, thisX, thisY, numRight[i] );
    thisX = xMin-1.0;
    if( numLeft[i] > 0.0 ) hs_fill_2d_hist( tmp_id, thisX, thisY, numLeft[i] );
  }

  delete [] numTop;
  delete [] numBottom;
  delete [] numRight;
  delete [] numLeft;

// Now set the overall statistics

// Histoscope doesn't want any of the sumW etc. We only need to
// worry about the number of entries. 

   hs_set_num_entries( tmp_id, numEntries );

  return temp;
}

void HepHistoHist2D::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. 

  if( temp.hid != 0 ) {
    char title[HS_MAX_TITLE_LENGTH];
    hs_title( histo_id_, title );
    int uid = 0;

// Histoscope's arithmetic functions for histogram objects will do the
// specified operation BUT puts the result in a new object and returns
// it's id. Find an unused category name to avoid clashes.

    char category[16];
    int * ids = new int[10];
    unsigned int numCategory = 1000; 
    int numMatches = 1;
    while( numMatches != 0 ) {
      sprintf(category,"%X",numCategory++);
      numMatches = hs_list_items( "", category, ids, 10, 1 );
    }
    delete [] ids;

    int hsid = hs_sum_histograms( uid, title, category,
                                  histo_id_, temp.hid, 1.0, 1.0 );
    hs_swap_data( hsid, histo_id_ );
    int numEntries = hs_num_entries( histo_id_ ) + hs_num_entries( temp.hid );
    hs_set_num_entries( histo_id_, numEntries );
    hs_delete( hsid ); 
  }
}

void HepHistoHist2D::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. 

  if( temp.hid != 0 ) {
    char title[HS_MAX_TITLE_LENGTH];
    hs_title( histo_id_, title );
    int uid = 0;

// Histoscope's arithmetic functions for histogram objects will do the
// specified operation BUT puts the result in a new object and returns
// it's id. Find an unused category name to avoid clashes.

    char category[16];
    int * ids = new int[10];
    unsigned int numCategory = 1000; 
    int numMatches = 1;
    while( numMatches != 0 ) {
      sprintf(category,"%X",numCategory++);
      numMatches = hs_list_items( "", category, ids, 10, 1 );
    }
    delete [] ids;

    int hsid = hs_sum_histograms( uid, title, category,
                                  histo_id_, temp.hid, 1.0, -1.0 );
    hs_swap_data( hsid, histo_id_ );
    int numEntries = hs_num_entries( histo_id_ ) + hs_num_entries( temp.hid );
    hs_set_num_entries( histo_id_, numEntries );
    hs_delete( hsid );
  }
}

void HepHistoHist2D::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. 

  if( temp.hid != 0 ) {
    char title[HS_MAX_TITLE_LENGTH];
    hs_title( histo_id_, title );
    int uid = 0;

// Histoscope's arithmetic functions for histogram objects will do the
// specified operation BUT puts the result in a new object and returns
// it's id. Find an unused category name to avoid clashes.

    char category[16];
    int * ids = new int[10];
    unsigned int numCategory = 1000; 
    int numMatches = 1;
    while( numMatches != 0 ) {
      sprintf(category,"%X",numCategory++);
      numMatches = hs_list_items( "", category, ids, 10, 1 );
    }
    delete [] ids;

    int hsid = hs_multiply_histograms( uid, title, category,
                                       histo_id_, temp.hid, 1.0 );
    hs_swap_data( hsid, histo_id_ );
    int numEntries = hs_num_entries( histo_id_ ) + hs_num_entries( temp.hid );
    hs_set_num_entries( histo_id_, numEntries );
    hs_delete( hsid );
  }
}

void HepHistoHist2D::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. 

  if( temp.hid != 0 ) {
    char title[HS_MAX_TITLE_LENGTH];
    hs_title( histo_id_, title );
    int uid = 0;

// Histoscope's arithmetic functions for histogram objects will do the
// specified operation BUT puts the result in a new object and returns
// it's id. Find an unused category name to avoid clashes.

    char category[16];
    int * ids = new int[10];
    unsigned int numCategory = 1000; 
    int numMatches = 1;
    while( numMatches != 0 ) {
      sprintf(category,"%X",numCategory++);
      numMatches = hs_list_items( "", category, ids, 10, 1 );
    }
    delete [] ids;

    int hsid = hs_divide_histograms( uid, title, category,
                                     histo_id_, temp.hid, 1.0 );
    hs_swap_data( hsid, histo_id_ );
    int numEntries = hs_num_entries( histo_id_ ) + hs_num_entries( temp.hid );
    hs_set_num_entries( histo_id_, numEntries );
    hs_delete( hsid );
  }
}

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

// If we got here, the temporary histogram object referred to by temp is
// known to be valid. 

  if( temp.hid != 0 ) hs_delete( temp.hid );

}

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

int   HepHistoHist2D::entries() const { // Number of entries in bins
  BADCALL(entries);

  // call histo function
  return hs_num_entries(histo_id_);
}

float HepHistoHist2D::weight() const {   // Sum of all weights
  BADCALL(weight);
  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_1d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);
  float sum = 0;
  int nBins = this->nBinsX()*this->nBinsY();
  for (int i=0;i<nBins;i++)
    sum += values_[i];
  return sum;
}

float HepHistoHist2D::weight2() const {   // Sum of all weights
  BADCALL(weight2);
  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_2d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);
  float sum = 0;
  int nBins = this->nBinsX()*this->nBinsY();
  for (int i=0;i<nBins;i++)
    sum += values_[i] * values_[i];
  return sum;
}

float HepHistoHist2D::sumX() const {  // retrieve weighted sum of X
  BADCALL(sumX);
  float sum = 0;
  float lowX = minX();
  int nBinsX = this->nBinsX();
  int nBinsY = this->nBinsY();
  float deltaX = (maxX()-lowX)/nBinsX;

  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_1d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);

  // loop over X
  for (int i=0;i<nBinsX;i++) {
    float theX = lowX + deltaX * (i+0.5);
    for (int j=0;j<nBinsY;j++)
      sum += theX * values_[ i*nBinsY + j ];
  } // end for i
  return sum;
}

float HepHistoHist2D::sumY() const {  // retrieve weighted sum of Y 
  BADCALL(sumY);
  float sum = 0;
  float lowY = minY();
  int nBinsX = this->nBinsX();
  int nBinsY = this->nBinsY();
  float deltaY = (maxY()-lowY)/nBinsY;

  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_1d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);

  // loop over X
  for (int i=0;i<nBinsY+0;i++) {
    float theY = lowY + deltaY * (i+0.5);
    for (int j=0;j<nBinsX;j++)
      sum += theY * values_[ j*nBinsX + i ];
  } // end for i
  return sum;
}

float HepHistoHist2D::sumX2() const {  // retrieve weighted sum of X^2
  BADCALL(sumX2);
  float sum = 0;
  float lowX = minX();
  int nBinsX = this->nBinsX();
  int nBinsY = this->nBinsY();
  float deltaX = (maxX()-lowX)/nBinsX;

  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_1d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);

  // loop over X
  for (int i=0;i<nBinsX;i++) {
    float theX = lowX + deltaX * (i+0.5);
    for (int j=0;j<nBinsY;j++)
      sum += theX * theX *  values_[ i*nBinsX + j ];
  } // end for i
  return sum;
}

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

  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_1d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);
 
  float* theY = new float[nBinsY+2];
  for (int k=0;k<nBinsY;k++)
    theY[k] = lowY + deltaY * (k+0.5);

  // loop over X
  for (int i=0;i<nBinsX;i++) {
    float theX = lowX + deltaX * (i+0.5);
    // loop over Y
    for (int j=0;j<nBinsY;j++) {
      sum += theY[j] * theX * values_[ i*nBinsX + j ];
    }
  } // end for i
  delete [] theY;
  return sum;

}

float HepHistoHist2D::sumY2() const {  // retrieve weighted sum of Y^2
  BADCALL(sumY2);
  float sum = 0;
  float lowY = minY();
  int nBinsX = this->nBinsX();
  int nBinsY = this->nBinsY();
  float deltaY = (maxY()-lowY)/nBinsY;

  // update our 'cached' values with the current content.
  // this copy is better than the only alternative: calling
  // hs_1d_hist_bin_value for each bin.
  hs_2d_hist_bin_contents(histo_id_, values_);
 
  // loop over Y
  for (int i=0;i<nBinsY;i++) {
    float theY = lowY + deltaY * (i+0.5);
    for (int j=0;j<nBinsX;j++)
      sum += theY * theY *values_[ j*nBinsX + i ];
  } // end for i
  return sum;
}

void HepHistoHist2D::reset()  {

  if ( mayUse() )  {
    hs_reset(histo_id_);
 }
  else  {
    ZMthrow(ZMxHepUnmanagedItem("attempt to reset() unmanaged histogram."));
    return;
  }
}

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