// ----------------------------------------------------------------------
//
// HepRootHistProf.cc - implementation of Root profile histogram
//
// HepRootHistProf is the implementation of the Root profile
// histogram. It is derived from HepHistProf (the generic interface) and
// HepRootHist. It is mainly a C++ wrapper around the  Root
// C++ histogramming libraries.
//
// No HepRootHistProf objects are ever directly
// instantiated by the programmer. When creating new HepHist* objects,
// the manager is in charge of creating the appropriate HepRoot object.
//
// To use generate Root histograms, HepRootFileManager must be used:
//
//   int lunit = 10;
//   HepFileManager* m = new HepRootFileManager("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);
//
//
//
//
//   12-Dec-1998  Philippe Canal  Creation
//   25-Feb-1999  J. Marraffino   Added increment operators for histograms
//
// ----------------------------------------------------------------------


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

#ifndef HEPROOTFILEMANAGER_H
#include "HepTuple/HepRootFileManager.h"
#endif

#ifndef HEPROOTHISTPROF_H
#include "HepTuple/HepRootHistProf.h"
#endif

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

#ifndef HEPTH1P_H
#include "HepTuple/HepTH1P.h"
#endif

#include "ZMutility/iostream"

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

#include <TH1.h>
#include <TProfile.h>

#include "HepTuple/HepTH1P.h"

ZM_BEGIN_NAMESPACE( zmht )	/*  namespace zmht  {  */

// Create a profile plot from an existing TProfile
HepRootHistProf::HepRootHistProf(
  HepRootFileManager * manager
  , TProfile* histo ) : 
  HepHistProf( manager, histo->GetTitle(), histo->GetNbinsX(),
	       histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax(),
// wrong for now:  histo->GetYaxis()->GetXmin(), histo->GetYaxis()->GetXmax(),
	       ((HepTH1P*)histo)->GetYmin(), ((HepTH1P*)histo)->GetYmax(),
	       histo->GetErrorOption(), histo->GetUniqueID() ),
  HepRootObj(histo)
{
  HEP_DEBUG( "HepRootHistProf::constructor( TProfile* " << histo << " )" );
}
 
// create a profile plot with the given title and binning
HepRootHistProf::HepRootHistProf(
  HepRootFileManager * 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.c_str(),
               nBinsX, lowX, highX, lowY, highY,
               chopt.c_str(), idReq )
{
  HEP_DEBUG( "HepRootHistProf::constructor( \"" << title << "\", " << *id_ << " )" );
  string name = mkname(title, *id_);

  // error= sigma/sqrt(N)
  data( new TProfile( name.c_str(), title.c_str(), nBinsX, 
		      lowX, highX, lowY, highY, chopt.c_str() ) );

  if ( data() == 0 ) {
    ZMthrow(ZMxHepCantMakeHist(string("Root failed to create the hist 1d: ") 
			       + title ) );
    isValid(false);
  }
}  // HepRootHistProf::HepRootHistProf()

HepRootHistProf::HepRootHistProf(
  HepRootFileManager * 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.c_str(), nBinsX, 
	       lowX, highX, lowY, highY, idReq )
{
  // Create a 1-D profile root histo with default error calculation:
  HEP_DEBUG( "HepRootHist1D::constructor( \"" << title << "\", " << id_ << " )" );
  
  string name = mkname(title, *id_);

  // error= sigma/sqrt(N)
  data( new TProfile( name.c_str(), title.c_str(), nBinsX, 
		      lowX, highX, lowY, highY ) );

  if ( data() == 0 ) {
    ZMthrow(ZMxHepCantMakeHist(string("Root failed to create the hist 1d: ") 
			       + title ) );
    isValid(false);
  } else {
     data()->SetUniqueID( idReq );
  }

}  // HepRootHistProf::HepRootHistProf()


HepRootHistProf::HepRootHistProf(
  HepRootHistProf * original
, const string& title
, const int     idReq
) : HepHistProf( original->manager(), title.c_str(),
                 original->nBinsX(), original->minX(), original->maxX(),
                 original->minY(), original->maxY(), idReq )
{

  data( (TNamed*)( original->data()->Clone() ) );
  if ( data() == 0 ) {
    ZMthrow(ZMxHepCantMakeHist(string("Root failed to create the hist 1d: ") 
			       + title ) );
    isValid(false);
  } else {
    data()->SetTitle( title.c_str() );
    data()->SetUniqueID( idReq );
  }

}  // HepRootHistProf::HepRootHistProf()


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


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

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

  if ( mayUse() )  {
    _accumulate( x, y, weight );
    ((TProfile*)data())->Fill( x, y, weight);
  }
  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;								\
  }

int   HepRootHistProf::entries() const { // number of entries
  CHECK(entries);

  // call histo function
  return (int) ((TH1*)data())->GetEntries();
}

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

float HepRootHistProf::weight() const {  // sum of all weight
  CHECK(weight);
  float result = 0.0;
  for(int i=1; i<nBinsX()+1; ++i) {
    result += weight(i);
  }
  return result;
}

float HepRootHistProf::weight2() const {  // sum of all weight
  CHECK(weight2);
  float result = 0.0;
  for(int i=1; i<nBinsX()+1; ++i) {
    result += weight(i)*weight(i);
  }
  return result;
}

float HepRootHistProf::weight( const int binNumX) const {  
  // sum of all weight in X-bin
  CHECK(weight); 
  return getWeight(binNumX);
}

float HepRootHistProf::sumX() const { // retrieve weighted sum of X
  if (! mayUse() ) {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to sum() an unmanaged histogram ")
	    +title()+string(".")));    
    return 0.0;
  }
  TH1 * hist = (TH1*)data();
  float result = 0.0;
  for(int i=1; i<nBinsX()+1; ++i) {
    result += (hist->GetBinCenter(i))*weight(i);
  }
  return result;
}

float HepRootHistProf::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;
  }
  TH1 * hist = (TH1*)data();
  float result = 0.0;
  float tempX;
  for(int i=1; i<nBinsX()+1; ++i) {
    tempX = hist->GetBinCenter(i);
    result += tempX*tempX*weight(i);
  }
  return result;
}
     
float HepRootHistProf::sumY() const {   // retrieve weighted sum of Y 
  CHECK(sumY);
  TProfile * hist = (TProfile*) data();
  float sumY = 0;
  for(int i=1; i<nBinsX()+1; ++i) {
    sumY += weight(i) * hist->GetBinContent(i);
  }
  return sumY;
}

float HepRootHistProf::sumY2() const {  // retrieve weighted sum of Y^2
  CHECK(sumY2);
  TProfile * hist = (TProfile*) data();
  float sumY2 = 0;
  for(int i=1; i<nBinsX()+1; ++i) {
    sumY2 += weight(i) * hist->GetBinContent(i) *
             weight(i) * hist->GetBinContent(i);
  }
  return sumY2;
}

float HepRootHistProf::sumY( const int binNumX) const { // retrieve sum of Y in X-bin
  CHECK(sumY);
  TProfile * hist = (TProfile*) data();
  if ( binNumX < 0 || binNumX > nBinsX() ) {
    return 0.0;
  }
  return  weight(binNumX) * hist->GetBinContent(binNumX);
}

float HepRootHistProf::sumY2( const int binNumX) const{ 
  // retrieve sum of Y^2 in X-bin
  CHECK(sumY2);
  TProfile * hist = (TProfile*) data();
  if ( binNumX < 0 || binNumX > nBinsX() ) {
    return 0.0;
  }
  return weight(binNumX) * hist->GetBinContent(binNumX) *
         weight(binNumX) * hist->GetBinContent(binNumX);
}

void HepRootHistProf::reset()  {

  if ( mayUse() )  {
    ((TH1*) data())->Reset();
 }
  else  {
    ZMthrow(ZMxHepUnmanagedItem("attempt to reset() unmanaged histogram."));
    return;
  }
}

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

  if ( mayUse() )  {
    return (float)(((TH1*)data())->GetBinContent(binNum));
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access to unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

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

  if ( mayUse() )  {
    return (float)(((TH1*)data())->GetBinError(binNum));
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access to unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

float HepRootHistProf::getWeight(	// get sum of weights for specified bin
  int numBin
) const  {

  if ( mayUse() )  {
    return ((HepTH1P*)data())->fetchWeight( numBin );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
    return 0.0;
  }
}

void HepRootHistProf::getSumW(	// get sum of weights per bin for all bins
  float * values
) const  {

  if ( mayUse() )  {
    ((HepTH1P*)data())->fetchSumW( values );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepRootHistProf::getSumWY( // get sum of weighted Y per bin for all bins
  float * values
) const  {

  if ( mayUse() )  {
    ((HepTH1P*)data())->fetchSumWY( values );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

void HepRootHistProf::getSumWY2( // get sum of weighted Y^2 per bin for all bins
  float * values
) const  {

  if ( mayUse() )  {
    ((HepTH1P*)data())->fetchSumWY2( values );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

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


void HepRootHistProf::getOverflows(		// get over/under-flow count
   float * numOver
 , float * numUnder
) const  {

  if ( mayUse() )  {
    *numUnder = this->bin( 0 );
    *numOver = this->bin( nBinsX()+1 );
  }
  else  {
    ZMthrow(ZMxHepUnmanagedItem(
	    string("attempt to access an unmanaged histogram ")
	    +title()+string(".")));
  }
}

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

// Manufacture a temporary Root histogram with all the characteristics of h.
// Make it a HepTH1P object so we can get to the protected data.

  HepTH1P * tmpHist = new HepTH1P( h );

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

  HepRawHist temp( tmpHist );
  return temp;

}

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

  TH1 * tmpHist = temp.HepTH1Pptr;
  if( tmpHist != 0 ) ((TH1*)data())->Add( tmpHist, 1.0 );

}

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

  TH1 * tmpHist = temp.HepTH1Pptr;
  if( tmpHist != 0 ) ((TH1*)data())->Add( tmpHist, -1.0 );

}

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

  ZMthrow(ZMxHepUnsupported("Root's Profile histograms have only a dummy Multiply method"));
  TH1 * tmpHist = temp.HepTH1Pptr;
  if( tmpHist != 0 ) ((TH1*)data())->Multiply( tmpHist );

}

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

  TH1 * tmpHist = temp.HepTH1Pptr;
  if( tmpHist != 0 ) ((TH1*)data())->Divide( tmpHist );

}

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

  HepTH1P * tmpHist = temp.HepTH1Pptr;
  if( tmpHist != 0 ) delete tmpHist;

}

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