// ----------------------------------------------------------------------
//
// HepHBookHistProf.cc - implementation of HBook profile histogram
//
// HepHBookHistProf is the implementation of the HBook profile
// histogram. It is derived from HepHistProf (the generic interface) and
// HepHBookHist. It is mainly a C++ wrapper around the CERNlib HBook
// Fortran histogramming libraries.
//
// No HepHBookHistProf 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);
//    :
//    :
//   HepHistProf* h(m->histProf("Title", 100, 0.0, 23.0, 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 HepHistProf. 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 HEPHBOOKHISTPROF_H
#include "HepTuple/HepHBookHistProf.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 profile plot with the given title and binning
HepHBookHistProf::HepHBookHistProf(
  HepHBookFileManager * manager
, const string& title
, const int nBinsX, const float lowX, const float highX
, const float lowY, const float highY
, const string& chopt
, const int idReq
) :
  HepHistProf( manager, title,
               nBinsX, lowX, highX, lowY, highY,
               chopt, idReq )
{
  HEP_DEBUG( "HepHBookHistProf::constructor( \"" << title << "\", " << *id_ << " )" );
}  // HepHBookHistProf::HepHBookHistProf()

HepHBookHistProf::HepHBookHistProf(
  HepHBookFileManager * manager
, const string& title
, const int nBinsX, const float lowX, const float highX
, const float lowY, const float highY
, const int idReq
) :
  HepHistProf( manager, title, nBinsX, lowX, highX, lowY, highY, idReq )
{
}  // HepHBookHistProf::HepHBookHistProf()


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

void HepHBookHistProf::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);
}

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


void HepHBookHistProf::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;
  }
}

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

float HepHBookHistProf::bin(                // retrieve contents of given bin
  const int binNum
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    int numBins = this->nBinsX();
    float result = 0.0;
    if( binNum > 0 && binNum <= numBins )  result = HI( *id_, binNum );
    HBookPopd();
    return result;
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access to unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

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

  if ( mayUse() )  {
    HBookPushd( dir() );
    int numBins = this->nBinsX();
    float result = 0.0;
    if( binNum > 0 && binNum <= numBins )  result = HIE( *id_, binNum );
    HBookPopd();
    return result;
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access to unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

void HepHBookHistProf::getSumW(		// get sum of weights for all bins
  float * data
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    char* type = "HIST";
    UNPAKW( *id_, *data, type, 1 );
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHistProf::getSumWY(	// get sum of weighted Y for all bins
  float * data
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    char* type = "HIST";
    UNPAKY( *id_, *data, type, 1 );
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepHBookHistProf::getSumWY2(	// get sum of weighted Y^2 for all bins
  float * data
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    char* type = "HIST";
    UNPAKYY( *id_, *data, type, 1 );
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

int   HepHBookHistProf::entries() const { // number of entries
  CHECK(entries);
  int result = 0;

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

  return result; 
}

int   HepHBookHistProf::entries( const int binNumX) const { 
  // number of entries in X-bin
  ZMthrow(ZMxHepUnsupported("HBook's Histogram Profile do not have access to the number of entries by bin"));
  CHECK(entries);
  return 0;
}

float HepHBookHistProf::weight() const {  // sum of all weight
  CHECK(weight);
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  float sumW = 0.0;

  HBookPushd( dir() );

  rtitle[maxTitleLen] = '\0';
  HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
  contents = new float[nx];
  char* type = "HIST";
  UNPAKW( *id_, *contents, type, 1 );
  for(int i=0; i<nx; ++i) {
    sumW += contents[i];
  }
  delete [] contents;
  HBookPopd();
  return sumW;
}

float HepHBookHistProf::weight2() const {  // sum of all weight
  CHECK(weight2);
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  float sumww = 0.0;

  HBookPushd( dir() );

  rtitle[maxTitleLen] = '\0';
  HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
  contents = new float[nx];
  char* type = "HIST";
  UNPAKW( *id_, *contents, type, 1 );
  HBookPopd();

  for(int i=0; i<nx; ++i) {
    sumww += contents[i] * contents[i];
  }
  delete [] contents;
  return sumww;
}

float HepHBookHistProf::weight( const int binNumX) const {  
  // sum of all weight in X-bin
  CHECK(weight); 
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  rtitle[maxTitleLen] = '\0';
  float sumW = 0.0;

  int numBins = this->nBinsX();
  if( binNumX > 0 && binNumX <= numBins ) {
    HBookPushd( dir() );
    HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
    contents = new float[nx];
    char* type = "HIST";
    UNPAKW( *id_, *contents, type, 1 );
    sumW = contents[binNumX-1];
    delete [] contents;
    HBookPopd();
  }
  return sumW;
}

float HepHBookHistProf::sumX() const { // retrieve weighted sum of X
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to sum() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  float sum = 0;
  float low = minX();
  int nBins = this->nBinsX();
  float delta = (maxX()-low)/nBins;
  HBookPushd( dir() );
  for (int i=1;i<nBins+1;i++)
    sum += (low + delta * (i-0.5)) * HI( *id_, i );
  HBookPopd();
  return sum;
}

float HepHBookHistProf::sumX2() const { // retrieve weighted sum of square (X^2)
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to sum() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  float sum = 0;
  float low = minX();
  int nBins = this->nBinsX();
  float delta = (maxX()-low)/nBins;
  HBookPushd( dir() );
  for (int i=1;i<nBins+1;i++) {
    float val = (low + delta * (i-0.5));
    sum += val * val * HI( *id_, i );
  }
  HBookPopd();
  return sum;
}
     
float HepHBookHistProf::sumY() const {   // retrieve weighted sum of Y 
  CHECK(sumY);
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  rtitle[maxTitleLen] = '\0';
  float sumY = 0.0;

  HBookPushd( dir() );
  HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
  contents = new float[nx];
  char* type = "HIST";
  UNPAKY( *id_, *contents, type, 1 );
  for(int i=0; i<nx; ++i) {
    sumY += contents[i];
  }
  delete [] contents;
  HBookPopd();
  return sumY;
}

float HepHBookHistProf::sumY2() const {  // retrieve weighted sum of Y^2
  CHECK(sumY2);
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  rtitle[maxTitleLen] = '\0';
  float sumY2 = 0.0;

  HBookPushd( dir() );
  HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
  contents = new float[nx];
  char* type = "HIST";
  UNPAKYY( *id_, *contents, type, 1 );
  for(int i=0; i<nx; ++i) {
    sumY2 += contents[i];
  }
  delete [] contents;
  HBookPopd();
  return sumY2;
}

float HepHBookHistProf::sumY( const int binNumX) const { // retrieve sum of Y in X-bin
  CHECK(sumY);
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  rtitle[maxTitleLen] = '\0';
  float sumY = 0.0;
//int kase = 2;

  int numBins = this->nBinsX();
  if( binNumX > 0 && binNumX <= numBins ) {
    HBookPushd( dir() );
    HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
    contents = new float[nx];
    char* type = "HIST";
    UNPAKY( *id_, *contents, type, 1 );
    sumY = contents[binNumX-1];
    delete [] contents;
    HBookPopd();
  }
  return sumY;
}

float HepHBookHistProf::sumY2( const int binNumX) const{ 
  // retrieve sum of Y^2 in X-bin
  CHECK(sumY2);
  int nx, ny, nwt, loc;
  float xmi, xma, ymi, yma;
  float *contents;
  const int maxTitleLen = 80;
  char rtitle[maxTitleLen+1];
  rtitle[maxTitleLen] = '\0';
  float sumY2 = 0.0;
//int kase = 3;

  int numBins = this->nBinsX();
  if( binNumX > 0 && binNumX <= numBins ) {
    HBookPushd( dir() );
    HGIVE( *id_, rtitle, nx, xmi, xma, ny, ymi, yma, nwt, loc );
    contents = new float[nx];
    char* type = "HIST";
    UNPAKYY( *id_, *contents, type, 1 );
    sumY2 = contents[binNumX-1];
    delete [] contents;
    HBookPopd();
  }
  return sumY2;
}

void HepHBookHistProf::reset()  {

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

void HepHBookHistProf::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 HepHBookHistProf::getOverflows(		// get over/under-flow count
   float * numOver
 , float * numUnder
) const  {

  if ( mayUse() )  {
    HBookPushd( dir() );
    int nBinsX = this->nBinsX();
    *numUnder = HI( *id_, 0 );
    *numOver = HI( *id_, nBinsX+1 );
    HBookPopd();
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

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

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

HepRawHist HepHBookHistProf::importHistProf( const HepHistProf & h ) {

// Gather up all the data we need from h

  int numBins = h.nBinsX();
  float minX = h.minX();
  float maxX = h.maxX();
  float * sumW   = new float[numBins];
  float * sumWY  = new float[numBins];
  float * sumWY2 = new float[numBins];

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

  h.getSumW  ( sumW );
  h.getSumWY ( sumWY );
  h.getSumWY2( sumWY2 );
  
  int numEntries;
  float sumWt;   
  float sumW2; 
  double sumWX;
  double sumWX2;
  h.getStatistics( numEntries, sumWt, 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 Prof histogram to park stuff

  char* title = "Temporary histProf for importHistProf";
  char* chopt = " ";
  HBPROF( id, title, numBins, h.minX(), h.maxY(),
                              h.minY(), h.maxY(), chopt );

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

  HepRawHist temp( id );

// Set up to move histogram data.

// Move the bin contents and bin errors

  char* type = "HIST";
  REPAKW ( id, *sumW,   type, 1 );
  REPAKY ( id, *sumWY,  type, 1 );
  REPAKYY( id, *sumWY2, type, 1 );
  delete [] sumW;
  delete [] sumWY;
  delete [] sumWY2;

// Then underflows and overflows

  if ( nover > 0 ) {
    float tempx = maxX+1.0;
    float tempy = 0.0;
    HFILL( id, tempx, tempy, nover );
  }
  if( nunder > 0 ) {
    float tempx = minX-1.0;
    float tempy = 0.0;
    HFILL( id, tempx, tempy, nunder );
  }

// Now set the overall statistics

  PUTSTAT( temp.hid, numEntries, sumWt, sumW2, sumWX, sumWX2 );

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

//HBookPopd();
  return temp;
}


void HepHBookHistProf::addHistProf( 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 importHistProf, 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 HepHBookHistProf::subtractHistProf( 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 importHistProf, 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 HepHBookHistProf::multiplyHistProf( 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 importHistProf, 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 HepHBookHistProf::divideHistProf( 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 importHistProf, 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 HepHBookHistProf::removeImportedHistProf( 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 importHistProf, 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  */