// ----------------------------------------------------------------------
//
// HepHistoHist1D.cc - implementation of Histo two-dimensional histogram
//
// HepHistoHist1D is the implementation of the Histo one-dimensional
// histogram. It is derived from HepHist1D (the generic interface) and
// HepHistoHist. It is mainly a C++ wrapper around the Histo
// C histogramming libraries.
//
// No HepHistoHist1D 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);
//    :
//    :
//   HepHist1D* h(m->hist1D("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 HepHist1D. 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, sum, sum2, entries
//   21-Feb-2000  Walter Brown   Improved C++ standard compliance
//
// ----------------------------------------------------------------------


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

#ifndef HEPHISTOHIST1D_H
#include "HepTuple/HepHistoHist1D.h"
#endif

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

#include "histoscope.h"

#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 1D histogram with the given title and binning
HepHistoHist1D::HepHistoHist1D(
  HepHistoFileManager * manager
, const string& title
, const int nBins, const float low, const float high
, const int idReq
) :
  HepHist1D( manager, title, nBins, low, high, idReq ), 
  values_( new float[nBins] )
{
  HEP_DEBUG( "HepHistoHist1D::constructor( \"" << title << "\", " << id_ << " )" );
  string name = title;
  histo_id_ = hs_create_1d_hist(*id_, (char*)title.c_str(), 
			  (char*)HepHistoFileManager::category(*dir_).c_str(),
			       (char*)(name+" x").c_str(),
			       (char*)(name+" y").c_str(),
			       nBins, low, high);
  if ( histo_id_ == 0 ) {
    ZMthrow(ZMxHepCantMakeHist(string("Histo failed to create the hist 1d: ") 
			       + title ) );
    isValid(false);
  }
}  //  HepHistoHist1D::HepHistoHist1D()

// create a 1D histogram with the given title and binning
HepHistoHist1D::HepHistoHist1D(
  HepHistoFileManager * manager
, const string& title
, const int nBins, const float low, const float high
, const int idReq
, const int histo_id
) :
  HepHist1D( manager, title, nBins, low, high, idReq ), 
  histo_id_(histo_id),
  values_( new float[nBins] )
{
  HEP_DEBUG( "HepHistoHist1D::constructor( \"" << title << "\", " << id_ << " )" );
}  //  HepHistoHist1D::HepHistoHist1D()

// clone an existing 1D histogram
HepHistoHist1D::HepHistoHist1D(
  HepHist1D * original
, const string&  title
, const int      idReq
) : HepHist1D( original->manager(), title,
               original->nBins(), original->min(), original->max(), idReq ),
    values_( new float[original->nBins()] )
{
  // make the Histoscope item
  int oldID = ((HepHistoHist1D*)original)->histo_id_;
  histo_id_ = hs_copy(oldID, idReq, (char*)title.c_str(), 
		 (char*)HepHistoFileManager::category(original->dir()).c_str() );

}  //  HepHistoHist1D::HepHistoHist1D()

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

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


void HepHistoHist1D::accumulate (
  const float x
, const float weight
)  {

  if ( mayUse() )  {
    _accumulate( x, weight );
    hs_fill_1d_hist( histo_id_, x, weight );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to accumulate() to unmanaged histogram ")
	    +title()+string(".")));
  }
}


float HepHistoHist1D::bin(                     // retrieve value of given bin
  const int binNum
) const  {

  if ( mayUse() )  {
    if ( binNum <= 0 || binNum > nBins() ) {
      float under,over;
      hs_1d_hist_overflows( histo_id_, &under, &over );
      if ( binNum <= 0 ) {
	return under;
      } else {
	return over;
      }
    } else {
      return hs_1d_hist_bin_value( histo_id_, binNum - 1 );
    }
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to bin() to unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

float HepHistoHist1D::binError(                // retrieve error of given bin
  const int binNum
) const  {

  if ( mayUse() )  {
    int numBins = nBins();
    float result;
    float * errors = new float[numBins];
    int status = hs_1d_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:
      ZMthrow(ZMxHepUnsupported(
			    string("binError: unknown error for histogram (")
			    +"1D"+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 HepHistoHist1D::getContents(             // get values for all bins
  float * data
) const  {

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

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

  if ( mayUse() )  {
    int numBins = nBins();
    float * errors = new float[numBins];
    int status = hs_1d_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:
      ZMthrow(ZMxHepUnsupported(
			    string("binError: unknown error for histogram (")
			    +"1D"+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
 *      get Statistics
 *	getErrors
 *	importHist1D
 *	addHist1D
 *	subtractHist1D
 *	multiplyHist1D
 *	divideHist1D
 *	removeImportedHist1D
 */

void HepHistoHist1D::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 numBins = nBins();
    float * contents = new float[numBins];
    getContents( contents );

    float thisBin;
    float thisX;
    float minX = min();
    float maxX = max();
    float binWidth = ( maxX - minX )/(float)numBins;
    for( int i = 0; i<numBins; ++i ) {
      thisBin = contents[i];
      thisX = minX + binWidth*((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 HepHistoHist1D::getOverflows(		// get over/under-flow count
   float * numOver
 , float * numUnder
) const  {

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

HepRawHist HepHistoHist1D::importHist1D( const HepHist1D & h ) {

// Gather up all the data we need from h

  int numBins = h.nBins();
  float minX = h.min();
  float maxX = h.max();
  float * contents  = new float[numBins];
  float * errors    = new float[numBins];

  float nunder;
  float nover;
  h.getOverflows( &nover, &nunder );

  h.getContents( contents );
  h.getErrors( errors );
  
  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 1D histogram to park stuff. Find an unused category
// name to avoid clashes.

  int uid = 0;
  int tmp_id = 0;
  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;

  tmp_id = hs_create_1d_hist( uid, "Temporary hist1D for importHist1D",
               category, "X", "Y", numBins, h.min(), h.max() );

// 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_1d_hist_block_fill( tmp_id, contents, errors, 0 );
  delete [] contents;
  delete [] errors;

// Then underflows and overflows

  if( nover  > 0 ) hs_fill_1d_hist( tmp_id, maxX+1.0,  nover );
  if( nunder > 0 ) hs_fill_1d_hist( tmp_id, minX-1.0, nunder );

// Now set the overall statistics

// Histoscope doesn't bother with 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 HepHistoHist1D::addHist1D( 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 always creates a new object to receive the sum histogram.
//  Let it do that and then put the id of the new object on top of the
//  one we have stored.  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 HepHistoHist1D::subtractHist1D( 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 always creates a new object to receive the sum histogram.
//  Let it do that and then put the id of the new object on top of the
//  one we have stored.  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 HepHistoHist1D::multiplyHist1D( 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 always creates a new object to receive the sum histogram.
//  Let it do that and then put the id of the new object on top of the
//  one we have stored.  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 HepHistoHist1D::divideHist1D( 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 always creates a new object to receive the sum histogram.
//  Let it do that and then put the id of the new object on top of the
//  one we have stored.  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 HepHistoHist1D::removeImportedHist1D( 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 HepHistoHist1D::entries() const { // number of entries
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to entries() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0;
  }
  // call histo function
  return hs_num_entries(histo_id_);
}

float HepHistoHist1D::weight() const { // sum of all weight
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to weight() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  
  // 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_1d_hist_bin_contents(histo_id_, values_);
  float sum = 0;
  int nBins = this->nBins();
  for (int i=0;i<nBins;i++)
    sum += values_[i];
  return sum;
}

float HepHistoHist1D::weight2() const { // sum of all weight
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to weight() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  
  // 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_1d_hist_bin_contents(histo_id_, values_);
  float sum = 0;
  int nBins = this->nBins();
  for (int i=0;i<nBins;i++)
    sum += values_[i] * values_[i];
  return sum;
}

float HepHistoHist1D::sum() const { // retrieve weighted sum
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to sum() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  // 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_1d_hist_bin_contents(histo_id_, values_);
  float sum = 0;
  float low = min();
  int nBins = this->nBins();
  float delta = (max()-low)/nBins;
  for (int i=0;i<nBins;i++)
    sum += (low + delta * (i+0.5)) * values_[i];
  return sum;
}

float HepHistoHist1D::sum2() const { // retrieve weighted sum of square (X^2)
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to sum() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  // 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_1d_hist_bin_contents(histo_id_, values_);
  float sum = 0;
  float low = min();
  int nBins = this->nBins();
  float delta = (max()-low)/nBins;
  for (int i=0;i<nBins;i++) {
    float val = (low + delta * (i+0.5));
    sum += val * val * values_[i];
  }
  return sum;
}


void HepHistoHist1D::reset()  {

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


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