#include "ATOOLS/Math/Histogram.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Org/MyStrStream.H"
#include "ATOOLS/Math/MathTools.H"
#include "ATOOLS/Org/My_File.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "ATOOLS/Org/My_MPI.H"
#include "ATOOLS/Org/CXXFLAGS.H"
#include "ATOOLS/Org/Scoped_Settings.H"
#include <stdio.h>

#include <sstream>

using namespace ATOOLS;

void Histogram::MPIInit()
{
#ifdef USING__MPI
  m_mfills=m_mpsfills=0.0;
  m_mvalues = new double*[m_depth];
  for (int j(0);j<m_depth;++j) {
    m_mvalues[j] = new double[m_nbin];
    for (int i(0);i<m_nbin;++i) m_mvalues[j][i]=0.0;
  }
#endif
}

Histogram::Histogram(int _type,double _lower,double _upper,int _nbin,
		     const std::string &name) :
  m_type(_type), m_nbin(_nbin), m_lower(_lower), m_upper(_upper), 
  m_yvalues(0),m_y2values(0), m_psvalues(0), m_tmp(0), m_fills(0), m_psfills(0), 
  m_finished(false), m_initialized(false), m_fuzzyexp(-1), m_name(name)
{
  m_ysums=NULL;
  m_ps2values=NULL;
  m_mcb = 0.;
  if (m_type>1000) {
    m_type-=1000;
    m_fuzzyexp = int(m_type/100);
    m_type-=100*m_fuzzyexp;
  }
  if (m_type>=100) {
    m_mcb = 1.;
    m_type-=100;
  }
  m_logarithmic = int(m_type/10);
  m_depth       = m_type-m_logarithmic*10+1;
  m_logbase = 1;
  switch(m_logarithmic) {
    case 1: 
      m_logbase = log(10.);
      m_upper   = log(m_upper)/m_logbase; m_lower = log(m_lower)/m_logbase;
      break;
    case 2: 
      m_upper   = log(m_upper); m_lower = log(m_lower);
      break;
    default: break;
  }
  m_binsize     = (m_upper-m_lower)/(double(m_nbin));

  if (m_binsize<=0.) {
    msg_Error()<<"Error in Histogram : Tried to initialize a histogram with  binsize <= 0 !"<<std::endl;
    m_active = 0;
    return;
  }

  m_nbin += 2;
  m_active = 1;
  m_yvalues   = new double[m_nbin];

  if (m_depth>1) {
    m_y2values   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_y2values[i]=0.;
    }
  }

  if (m_depth>2) {
    m_psvalues   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_psvalues[i]=0.;
    }
  }
  if (m_depth>3) {
    m_ps2values   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_ps2values[i]=0.;
    }
  }

  if (m_mcb!=0.) {
    m_tmp   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_tmp[i]=0.;
    }
  }

  for (int i=0;i<m_nbin;i++) {
    m_yvalues[i]=0.;
  }
  MPIInit();
}

Histogram::Histogram(const Histogram * histo)
  : m_yvalues(0), m_y2values(0), m_psvalues(0), m_tmp(0) {
  m_ysums=NULL;
  CopyFrom(histo);
}

void Histogram::CopyFrom(const Histogram *histo)
{
  if (m_yvalues) delete [] m_yvalues;
  if (m_y2values) delete [] m_y2values;
  if (m_psvalues) delete [] m_psvalues;
  if (m_ysums) delete [] m_ysums;
  if (m_tmp) delete [] m_tmp;

  m_lower   = histo->m_lower;
  m_upper   = histo->m_upper;
  m_logbase = histo->m_logbase;
  m_logarithmic = histo->m_logarithmic;
  m_mcb     = histo->m_mcb;
  m_nbin    = histo->m_nbin;
  m_depth   = histo->m_depth;
  m_type    = histo->m_type;
  m_fills   = histo->m_fills;
  m_psfills = histo->m_psfills;

  m_binsize = histo->m_binsize;
  m_active  = 1;
  m_finished = histo->m_finished;
  m_name     = histo->m_name;

  m_yvalues   = new double[m_nbin];
  for (int i=0;i<m_nbin;i++) {
    m_yvalues[i]  = histo->m_yvalues[i]; 
  }
  if (m_depth>1) {
    m_y2values   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_y2values[i]=histo->m_y2values[i];
    }
  }
  if (m_depth>2) {
    m_psvalues   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_psvalues[i]=histo->m_psvalues[i];
    }
  }
  if (m_mcb!=0.) {
    m_tmp   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_tmp[i]=0.;
    }
  }
  if (histo->m_ysums==NULL) m_ysums=NULL;
  else {
    m_ysums = new double[m_nbin];
    for (int i=0;i<m_nbin;i++)
      m_ysums[i]=histo->m_ysums[i];
  }
  MPIInit();
}

Histogram::Histogram(const std::string & pID,const int mode,std::string content)
  :  m_yvalues(0), m_y2values(0), m_psvalues(0), m_tmp(0), m_fills(0), m_mcb(0.)  {
  m_ysums=NULL;
  m_finished=true;

  std::stringstream ifile;
 
  if (content=="") {
    My_In_File file(pID.c_str());
    file.Open();
    std::istream *is=dynamic_cast<std::istream*>(&*file);
    ifile << (*is).rdbuf();
    file.Close();    
  }
  else {
    ifile.str(content);
  }

    ifile>>m_type>>m_nbin>>m_lower>>m_upper;


    m_logarithmic = int(m_type/10);
    m_depth       = m_type-m_logarithmic*10+1;

    m_logbase = 1;
    switch(m_logarithmic) {
    case 1: 
      m_logbase = log(10.);
      break;
    default: break;
    }
    m_binsize     = (m_upper-m_lower)/(double(m_nbin-2));
    
    if (m_binsize<=0.) {
      msg_Error()<<"Error in Histogram : "
		 <<"Tried to initialize a histogram with m_binsize <= 0 !"
		 <<std::endl;
      m_active = 0;
      return;
    }
    
    std::string helps;
    m_active = 1;
    m_yvalues   = new double[m_nbin];
    if (mode==0) ifile>>m_yvalues[0];
    else {
      ifile>>helps;
      if (helps.find("nan")!=std::string::npos) helps="0";
      m_yvalues[0]=ToType<double>(helps);
    }
    if (m_depth>1) {
      m_y2values   = new double[m_nbin];
      if (mode==0) ifile>>m_y2values[0];
      else {
	ifile>>helps;
	if (helps.find("nan")!=std::string::npos) helps="0";
	m_y2values[0]=ToType<double>(helps);
      }
    }    
    if (m_depth>2) {
      m_psvalues   = new double[m_nbin];
      if (mode==0) ifile>>m_psvalues[0];
      else {
	ifile>>helps;
	if (helps.find("nan")!=std::string::npos) helps="0";
	m_psvalues[0]=ToType<double>(helps);
      }
    }    

    if (mode==0) ifile>>m_yvalues[m_nbin-1];
    else {
      ifile>>helps;
      if (helps.find("nan")!=std::string::npos) helps="0";
      m_yvalues[m_nbin-1]=ToType<double>(helps);
    }
    if (m_depth>1) {
      if (mode==0) ifile>>m_y2values[m_nbin-1];
      else {
	ifile>>helps;
	if (helps.find("nan")!=std::string::npos) helps="0";
	m_y2values[m_nbin-1]=ToType<double>(helps);
      }
    }    
    if (m_depth>2) {
      if (mode==0) ifile>>m_psvalues[m_nbin-1];
      else {
	ifile>>helps;
	if (helps.find("nan")!=std::string::npos) helps="0";
	m_psvalues[m_nbin-1]=ToType<double>(helps);
      }
    }    
    ifile>>m_fills;

  double x;
  for (int i=0;i<m_nbin-1;i++) {
    ifile>>x;
    if (!IsEqual(x,m_lower+i*m_binsize,ifile.precision()-1)) {
      msg_Error()<<METHOD<<"(): Corrupted input file '"<<pID<<"'."<<std::endl;
      m_active=0;
      break;
    }
    if (mode==0) ifile>>m_yvalues[i+1];
    else {
      ifile>>helps;
      if (helps.find("nan")!=std::string::npos) helps="0";
      m_yvalues[i+1]=ToType<double>(helps);
    }
    if (m_depth>1) {
      if (mode==0) ifile>>m_y2values[i+1];
      else {
	ifile>>helps;
	if (helps.find("nan")!=std::string::npos) helps="0";
	m_y2values[i+1]=ToType<double>(helps);
      }
      m_y2values[i+1] = sqr(m_y2values[i+1]);
    }    
    if (m_depth>2) {
      if (mode==0) ifile>>m_psvalues[i+1];
      else {
	ifile>>helps;
	if (helps.find("nan")!=std::string::npos) helps="0";
	m_psvalues[i+1]=ToType<double>(helps);
      }
    }    
  }
  if (ifile.eof()) {
    msg_Error()<<METHOD<<"(): Corrupted input file '"<<pID<<"'."<<std::endl;
    m_active=0;
  }
  ifile>>x;
  if (!ifile.eof()) {
    msg_Error()<<METHOD<<"(): Corrupted input file '"<<pID<<"'."<<std::endl;
    m_active=0;
  }
  MPIInit();
}


Histogram::~Histogram() {
    
  if (m_yvalues!=0) { 
    delete [] m_yvalues; m_yvalues = 0; 
  }
  if (m_y2values!=0) { 
    delete [] m_y2values; m_y2values = 0; 
  }
  if (m_psvalues!=0) { 
    delete [] m_psvalues; m_psvalues = 0; 
  }
  if (m_tmp!=0) { 
    delete [] m_tmp; m_tmp = 0; 
  }
  if (m_ysums!=0) delete [] m_ysums;
#ifdef USING__MPI
  for (int j(0);j<m_depth;++j) delete [] m_mvalues[j];
  delete [] m_mvalues;
#endif
}


void Histogram::Finalize() {
  if (!m_finished) {
    m_finished=true;
    if (m_fills==0.) return;
    for (int i=0;i<m_nbin;++i) {
      m_yvalues[i]/=m_fills*m_binsize;
      if (m_depth>1) {
 	m_y2values[i]/=m_fills*sqr(m_binsize);
	if (m_fills>1) m_y2values[i]=(m_y2values[i]-sqr(m_yvalues[i]))/(m_fills-1);
      }
    }
    if (m_depth>2) {
      double itg = Integral()/(m_psfills*m_binsize);
      for (int i=0;i<m_nbin;++i) {
	m_psvalues[i]*=itg;
      }
    }
  }
}

void Histogram::Restore() {
  if (m_finished) {
    for (int i=0;i<m_nbin;++i) {
      if (m_depth>1) {
	if (m_fills>1) m_y2values[i]=(m_fills-1)*m_y2values[i]+sqr(m_yvalues[i]);
	m_y2values[i]*=m_fills*sqr(m_binsize);
	if (m_depth>2) {
	  m_psvalues[i]*=m_psfills*m_binsize;
	}
      }
      m_yvalues[i]*=m_fills*m_binsize;
    }
    m_finished=false;
  }
}

double Histogram::Average() const
{
  double prod=0., sum=0., bin=m_lower;
  double width=(m_upper-m_lower)/m_nbin;
  bin += width/2.;

  for (int i=1;i<m_nbin-1;++i) {
    prod += m_yvalues[i]*bin;
    bin  += width;
    sum  += m_yvalues[i];
  }
  return prod/sum;
}

double Histogram::Mean() const
{
  double sum=0., range=0.;
  for (int i=1;i<m_nbin-1;++i) {
    double width=(m_upper-m_lower)/m_nbin;
    if (m_logarithmic) 
      width=pow(m_logbase,m_lower+i*width)-pow(m_logbase,m_lower+(i-1)*width);  
    sum+=m_yvalues[i]*width;
    range+=width;
  }
  return sum/range;
}

void Histogram::Reset() {
  for (int i=0;i<m_nbin;i++) { 
    m_yvalues[i]=0.;
    if (m_depth>1) {
      m_y2values[i]=0.;
    }
    if (m_depth>2) {
      m_psvalues[i]=0.;
    }
  }
  m_fills=0;
  m_psfills=0;
}

void Histogram::Scale(double scale) {
  for (int i=0;i<m_nbin;i++) { 
    m_yvalues[i]*= scale;
    if (m_depth>1) m_y2values[i]*=sqr(scale); 
    if (m_depth>2) m_psvalues[i]*=scale; 
  }
}

void Histogram::ScaleHistogramWidth(double factor, int mode)
{
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a "
	       <<"histogram with binsize <= 0 !"<<std::endl;
    return;
  }
  if (factor<=0) {
    msg_Error()<<"Error in Histogram : Tried to scale binsize "
               <<"of a histogram with a factor <= 0 !"<<std::endl;
    return;
  }
  m_upper   *= factor;
  m_binsize *= factor;
  if (!mode) {
    for (int i=0;i<m_nbin;i++) {
      m_yvalues[i] /= factor;
      if (m_depth>1) m_y2values[i] /= sqr(factor); 
      if (m_depth>2) m_psvalues[i] /= factor;
    }
  }
}

void Histogram::Output() {
  if (!msg_LevelIsDebugging()) return;
  msg_Out()<<"----------------------------------------"<<std::endl
	   <<"    "<<m_yvalues[0]<<std::endl
	   <<"----------------------------------------"<<std::endl;
  double result = 0.;
  for (int i=0;i<m_nbin-2;i++) {
    msg_Out()<<m_lower+i*m_binsize<<"  ";
    msg_Out()<<m_yvalues[i+1]<<"  ";
    if (m_depth>1) msg_Out()<<sqrt(m_y2values[i+1]);
    result += m_yvalues[i+1];
    msg_Out()<<std::endl;
  }
  msg_Out()<<m_lower+(m_nbin-2)*m_binsize<<" == "<< m_upper<<std::endl
	   <<"----------------------------------------"<<std::endl
	   <<"    "<<m_yvalues[m_nbin-1]<<std::endl
	   <<"----------------------------------------"<<std::endl
	   <<"Inside the range : "<<result<<std::endl;
}

double Histogram::LowEdge(int i) const 
{
  if (m_logarithmic) return exp(m_logbase*(m_lower+i*m_binsize));  
  return m_lower+i*m_binsize;
}

double Histogram::HighEdge(int i) const 
{
  if (m_logarithmic) return exp(m_logbase*(m_lower+(i+1)*m_binsize));  
  return m_lower+(i+1)*m_binsize;
}


void Histogram::Output(const std::string name) 
{
  if (!m_active) return;
  My_Out_File ofile(name);
  ofile.Open();
  if (rpa) {
    Settings& s = Settings::GetMainSettings();
    ofile->precision(
        s["HISTOGRAM_OUTPUT_PRECISION"].SetDefault(6).GetScalar<int>());
  }

  if (m_fills>=0) {
    *ofile<<m_type<<" "<<m_nbin<<" "<<m_lower<<" "<<m_upper<<" ";
    *ofile<<m_yvalues[0]<<"  ";
    if (m_depth>1) *ofile<<m_y2values[0]<<"  ";
    *ofile<<m_yvalues[m_nbin-1]<<"  ";
    if (m_depth>1) *ofile<<m_y2values[m_nbin-1]<<"  ";
    *ofile<<m_fills<<"\n";
  }
  for (int i=0;i<m_nbin-1;i++) {
    *ofile<<m_lower+i*m_binsize<<"  ";
    *ofile<<m_yvalues[i+1]<<"  ";
    if (m_depth>1) *ofile<<sqrt(m_y2values[i+1])<<"  ";
    if (m_depth>2) *ofile<<m_psvalues[i+1]<<"  ";
    if (m_depth>3) *ofile<<sqrt(m_ps2values[i+1])<<"  ";
    *ofile<<"\n";
  }
  ofile.Close();
}

void Histogram::MPISync()
{
#ifdef USING__MPI
  int size=mpi->Size();
  if (size>1) {
    int cn=m_depth*m_nbin+2;
    double *values = new double[cn];
    for (int j(0);j<m_depth;++j)
      for (int i(0);i<m_nbin;++i) values[j*m_nbin+i]=m_mvalues[j][i];
    values[cn-2]=m_mfills;
    values[cn-1]=m_mpsfills;
    mpi->Allreduce(values,cn,MPI_DOUBLE,MPI_SUM);
    for (int j(0);j<m_depth;++j)
      for (int i(0);i<m_nbin;++i) m_mvalues[j][i]=values[j*m_nbin+i];
    m_mfills=values[cn-2];
    m_mpsfills=values[cn-1];
    delete [] values;
  }
  for (int i(0);i<m_nbin;++i) {
    m_yvalues[i]+=m_mvalues[0][i];
    m_mvalues[0][i]=0.0;
    if (m_depth>1) {
      m_y2values[i]+=m_mvalues[1][i];
      m_mvalues[1][i]=0.0;
    }
    if (m_depth>2) {
      m_psvalues[i]+=m_mvalues[2][i];
      m_mvalues[2][i]=0.0;
    }
    if (m_depth>3) {
      m_ps2values[i]+=m_mvalues[3][i];
      m_mvalues[3][i]=0.0;
    }
  }
  m_fills+=m_mfills;
  m_psfills+=m_mpsfills;
  m_mfills=m_mpsfills=0.0;
#endif
}

double Histogram::GeneratePoint(const double &rn)
{
  if (m_ysums==NULL) {
    m_ysums = new double[m_nbin];
    double sum=0.;
    for (int i=0;i<m_nbin;i++)
      m_ysums[i]=sum+=m_yvalues[i];
  }
  double disc=rn*m_ysums[m_nbin-1];
  int l(0), r(m_nbin-1), c((l+r)/2);
  double a(m_ysums[c]);
  while (r-l>1) {
    if (disc<a) r=c;
    else l=c;
    c=(l+r)/2;
    a=m_ysums[c];
  }
  double x(0.);
  if (disc<m_ysums[l]) x=m_lower+(l-1+disc/m_yvalues[l])*m_binsize;
  else x=m_lower+(r-1+(disc-m_ysums[l])/m_yvalues[r])*m_binsize;
  if (m_logarithmic>0) x=exp(m_logbase*x);
  return x;
}

void Histogram::Insert(double coordinate) {
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a "
	       <<"histogram with binsize <= 0 !"<<std::endl;
    return;
  }

#ifdef USING__MPI
  m_mfills++;

  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;
  if (coordinate<m_lower) { m_mvalues[0][0]       += double(1); return; }
  if (coordinate>m_upper) { m_mvalues[0][m_nbin-1] += double(1); return; }
  for (int i=1;i<m_nbin-1;i++) {
    if ( (coordinate >= m_lower + (i-1)*m_binsize) &&
	 (coordinate <  m_lower + i*m_binsize) ) {
      m_mvalues[0][i] += double(1); 
      return; 
    }
  }
#else
  m_fills++;

  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;
  if (coordinate<m_lower) { m_yvalues[0 ]       += double(1); return; }
  if (coordinate>m_upper) { m_yvalues[m_nbin-1] += double(1); return; }
  for (int i=1;i<m_nbin-1;i++) {
    if ( (coordinate >= m_lower + (i-1)*m_binsize) &&
	 (coordinate <  m_lower + i*m_binsize) ) {
      m_yvalues[i] += double(1); 
      return; 
    }
  }
#endif
}

void Histogram::Insert(int i,double value,double ncount) {
  if (IsBad(value)) {
    msg_Error()<<METHOD<<"("<<i<<","<<value<<","<<ncount<<"): Skip bad weight.";
    return;
  }
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a "
			  <<"histogram with binsize <= 0 !"<<std::endl;
    return;
  }
#ifdef USING__MPI
  m_mfills+=ncount;
  if (value==0.) return;
  m_mpsfills++;

  if (i<0) { 
    m_mvalues[0][0] += value;
    if (m_depth>1) {
      if (value>m_mvalues[1][0]) m_mvalues[1][0] = value;
      if (m_depth>2) m_mvalues[2][0] += 1.;
    }
    return; 
  }

  if (i>=m_nbin) { 
    m_mvalues[0][m_nbin-1] += value; 
    if (m_depth>1) {
      if (value>m_mvalues[1][m_nbin-1]) m_mvalues[1][m_nbin-1] = value;
      if (m_depth>2) m_mvalues[2][m_nbin-1] += 1.;
    }
    return; 
  }

  m_mvalues[0][i] += value;
  if (m_depth>1) {
    m_mvalues[1][i] += value*value;
    if (m_depth>2) m_mvalues[2][i] += 1.;
  }
#else
  m_fills+=ncount;
  if (value==0.) return;
  m_psfills++;

  if (i<0) { 
    m_yvalues[0] += value;
    if (m_depth>1) {
      if (value>m_y2values[0]) m_y2values[0] = value;
      if (m_depth>2) m_psvalues[0] += 1.;
    }
    return; 
  }

  if (i>=m_nbin) { 
    m_yvalues[m_nbin-1] += value; 
    if (m_depth>1) {
      if (value>m_y2values[m_nbin-1]) m_y2values[m_nbin-1] = value;
      if (m_depth>2) m_psvalues[m_nbin-1] += 1.;
    }
    return; 
  }

  m_yvalues[i] += value;
  if (m_depth>1) {
    m_y2values[i] += value*value;
    if (m_depth>2) m_psvalues[i] += 1.;
  }
#endif
}

void Histogram::InsertMCB(double coordinate,double value,double ncount) {
  if (IsBad(value)) {
    msg_Error()<<METHOD<<"("<<coordinate<<","<<value<<","<<ncount<<"): Skip bad weight.";
    return;
  }
  if (!m_tmp) {
    m_tmp   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_tmp[i]=0.;
    }
  }
  m_mcb = ncount;

  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

  int bin = int((coordinate-m_lower)/m_binsize+1.);
  if (bin<0) bin=0;
  if (bin>=m_nbin) bin=m_nbin-1;
  if (bin==0||bin==m_nbin-1) {
    m_tmp[bin] += value;
    return;
  }

  double x = (coordinate-m_lower)/m_binsize-double(bin)+0.5;
  if ((bin==1&&x<0.)||(bin==m_nbin-2&&x>0.)) {
    m_tmp[bin] += value;
    return;
  }
  double ff=1.;
  if (m_fuzzyexp==0) ff=0.5;
  if (m_fuzzyexp>0) ff=1.-0.5*pow(2.*dabs(x),m_fuzzyexp);
  if (m_fuzzyexp==9) ff=1.-0.5*sqrt(2.*dabs(x));

  m_tmp[bin] += ff*value; 
  if (x>0.) m_tmp[bin+1] += (1.-ff)*value;
  if (x<0.) m_tmp[bin-1] += (1.-ff)*value;
}

void Histogram::InsertMCBIM(double coordinate,double value) {
  if (IsBad(value)) {
    msg_Error()<<METHOD<<"("<<coordinate<<","<<value<<"): Skip bad weight.";
    return;
  }
  if (!m_tmp) {
    m_tmp   = new double[m_nbin];
    for (int i=0;i<m_nbin;i++) {
      m_tmp[i]=0.;
    }
  }
  m_mcb = 1.;
  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

  int bin = int((coordinate-m_lower)/m_binsize+1.);
  if (bin<0) bin=0;
  if (bin>=m_nbin) bin=m_nbin-1;
  for (int i=bin+1;i<m_nbin;i++) m_tmp[i] += value;
}

void Histogram::FinishMCB()
{
  if (!m_tmp) {
    m_tmp = new double[m_nbin];
    for (int i=0;i<m_nbin;++i) m_tmp[i]=0.0;
  }
#ifdef USING__MPI
  m_mfills+=m_mcb;
  m_mpsfills++;

  for (int i=0;i<m_nbin;i++) {
    m_mvalues[0][i] += m_tmp[i];
    if (m_depth>1) {
      m_mvalues[1][i] += m_tmp[i]*m_tmp[i];
      if (m_depth>2&&m_tmp[i]!=0.) m_mvalues[2][i] += 1.;
    }
    m_tmp[i] = 0.;
  }
#else
  m_fills+=m_mcb;
  m_psfills++;

  for (int i=0;i<m_nbin;i++) {
    m_yvalues[i] += m_tmp[i];
    if (m_depth>1) {
      m_y2values[i] += m_tmp[i]*m_tmp[i];
      if (m_depth>2&&m_tmp[i]!=0.) m_psvalues[i] += 1.;
    }
    m_tmp[i] = 0.;
  }
#endif
}

void Histogram::Insert(double coordinate,double value,double ncount) {
  if (IsBad(value)) {
    msg_Error()<<METHOD<<"("<<coordinate<<","<<value<<","<<ncount<<"): Skip bad weight.";
    return;
  }
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a "
			  <<"histogram with binsize <= 0 !"<<std::endl;
    return;
  }
#ifdef USING__MPI
  m_mfills+=ncount;
  if (value==0.) return;
  m_mpsfills++;

  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

  int bin = int((coordinate-m_lower)/m_binsize+1.);
  if (bin<0) bin=0;
  if (bin>=m_nbin) bin=m_nbin-1;
  if (bin==0||bin==m_nbin-1) {
    m_mvalues[0][bin] += value;
    if (m_depth>1) {
      if (value>m_mvalues[1][bin]) m_mvalues[1][bin] = value;
      if (m_depth>2) m_mvalues[2][bin] += 1.;
    }
    return;
  }

  m_mvalues[0][bin] += value; 
  if (m_depth>1) {
    m_mvalues[1][bin] += sqr(value);
    if (m_depth>2) m_mvalues[2][bin] += 1.;
  }
 
  if (m_fuzzyexp<0) return;

  double x = (coordinate-m_lower)/m_binsize-double(bin)+0.5;
  if (bin==1&&x<0.) return;
  if (bin==m_nbin-2&&x>0.) return;
  double ff=1.;
  if (m_fuzzyexp==0) ff=0.5;
  if (m_fuzzyexp>0) ff=0.5*pow(2.*dabs(x),m_fuzzyexp);
  if (m_fuzzyexp==9) ff=0.5*sqrt(2.*dabs(x));

  m_mvalues[0][bin] -= ff*value; 
  if (m_depth>1) {
    m_mvalues[1][bin] += sqr(ff*(value))-sqr(value);
    if (m_depth>2) m_mvalues[2][bin] -= ff;
  }

  if (x>0.) {
    m_mvalues[0][bin+1] += ff*value;
    if (m_depth>1) {
      m_mvalues[1][bin+1] += sqr(ff*value);
      if (m_depth>2) m_mvalues[2][bin+1] += ff;
    }
  }
  if (x<0.) {
    m_mvalues[0][bin-1] += ff*value;
    if (m_depth>1) {
      m_mvalues[1][bin-1] += sqr(ff*value);
      if (m_depth>2) m_mvalues[2][bin-1] += ff;
    }
  }
#else
  m_fills+=ncount;
  if (value==0.) return;
  m_psfills++;

  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

  int bin = int((coordinate-m_lower)/m_binsize+1.);
  if (bin<0) bin=0;
  if (bin>=m_nbin) bin=m_nbin-1;
  if (bin==0||bin==m_nbin-1) {
    m_yvalues[bin] += value;
    if (m_depth>1) {
      if (value>m_y2values[bin]) m_y2values[bin] = value;
      if (m_depth>2) m_psvalues[bin] += 1.;
    }
    return;
  }

  m_yvalues[bin] += value; 
  if (m_depth>1) {
    m_y2values[bin] += sqr(value);
    if (m_depth>2) m_psvalues[bin] += 1.;
  }
 
  if (m_fuzzyexp<0) return;

  double x = (coordinate-m_lower)/m_binsize-double(bin)+0.5;
  if (bin==1&&x<0.) return;
  if (bin==m_nbin-2&&x>0.) return;
  double ff=1.;
  if (m_fuzzyexp==0) ff=0.5;
  if (m_fuzzyexp>0) ff=0.5*pow(2.*dabs(x),m_fuzzyexp);
  if (m_fuzzyexp==9) ff=0.5*sqrt(2.*dabs(x));

  m_yvalues[bin] -= ff*value; 
  if (m_depth>1) {
    m_y2values[bin] += sqr(ff*(value))-sqr(value);
    if (m_depth>2) m_psvalues[bin] -= ff;
  }

  if (x>0.) {
    m_yvalues[bin+1] += ff*value;
    if (m_depth>1) {
      m_y2values[bin+1] += sqr(ff*value);
      if (m_depth>2) m_psvalues[bin+1] += ff;
    }
  }
  if (x<0.) {
    m_yvalues[bin-1] += ff*value;
    if (m_depth>1) {
      m_y2values[bin-1] += sqr(ff*value);
      if (m_depth>2) m_psvalues[bin-1] += ff;
    }
  }
#endif
}

void Histogram::InsertRange(double start, double end, double value) {
  if (IsBad(value)) {
    msg_Error()<<METHOD<<"("<<start<<","<<end<<","<<value<<"): Skip bad weight.";
    return;
  }
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a "
			  <<"histogram with binsize <= 0 !"<<std::endl;
    return;
  }
  if (m_logarithmic>0) {
    if (start>0)
      start = log(start)/m_logbase;
    else
      start = -30;
    if (end>0)
      end = log(end)/m_logbase;
    else 
      end = -30;
  }
#ifdef USING__MPI
  m_mfills++;

  // underrun
  if (start<m_lower) { 
    m_mvalues[0][0] += value;
    if (m_depth>1) {
      if (value>m_mvalues[1][0]) m_mvalues[1][0] = value;
    }
  }

  // overflow
  if (start>m_upper) { 
    m_mvalues[0][m_nbin-1] += value; 
    if (m_depth>1) {
      if (value>m_mvalues[1][m_nbin-1]) m_mvalues[1][m_nbin-1] = value;
    }
  }

  double low,up;
  int hit=0;
  low = m_lower; up = m_lower+m_binsize;
  for (int i=1;i<m_nbin-1;i++) {
    if ((start < up) && (end >= low) ) {
      double fac=1;
      if ((start<=low)&&(up<=end )) {
	m_mvalues[0][i] += value;
      } 
      else if ((low<start)&&(up<=end)) {
	fac = (up-start)/m_binsize;
	m_mvalues[0][i] += value *fac;
      }
      else if ((start<=low)&&(end < up)) {
	fac = (end-low)/m_binsize;
	m_mvalues[0][i] += value *fac;
      }
      else if ((low<start)&&(end <up)) {
	fac = (end-start)/m_binsize;
	m_mvalues[0][i] += value *fac;
      }
    

      hit=1;
      if (m_depth>1) {
	if (value>m_mvalues[1][i]) m_mvalues[1][i] = value;
      }
    }
    low = up;
    up += m_binsize;
  }
#else
  m_fills++;

  // underrun
  if (start<m_lower) { 
    m_yvalues[0] += value;
    if (m_depth>1) {
      if (value>m_y2values[0]) m_y2values[0] = value;
    }
  }

  // overflow
  if (start>m_upper) { 
    m_yvalues[m_nbin-1] += value; 
    if (m_depth>1) {
      if (value>m_y2values[m_nbin-1]) m_y2values[m_nbin-1] = value;
    }
  }

  double low,up;
  low = m_lower; up = m_lower+m_binsize;
  for (int i=1;i<m_nbin-1;i++) {
    if ((start < up) && (end >= low) ) {
      double fac=1;
      if ((start<=low)&&(up<=end )) {
	m_yvalues[i] += value;
      } 
      else if ((low<start)&&(up<=end)) {
	fac = (up-start)/m_binsize;
	m_yvalues[i] += value *fac;
      }
      else if ((start<=low)&&(end < up)) {
	fac = (end-low)/m_binsize;
	m_yvalues[i] += value *fac;
      }
      else if ((low<start)&&(end <up)) {
	fac = (end-start)/m_binsize;
	m_yvalues[i] += value *fac;
      }
    
      if (m_depth>1) {
	if (value>m_y2values[i]) m_y2values[i] = value;
      }
    }
    low = up;
    up += m_binsize;
  }
#endif
}



double Histogram::Bin(int bin) const { return m_yvalues[bin]; }
double Histogram::Bin2(int bin) const { return m_y2values[bin]; }

double Histogram::Bin(double coordinate) const
{ 
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a histogram wit binsize <= 0 ! Return 0.."<<std::endl;
    return -1.0;
  }
  else {
    if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

    if (coordinate<m_lower) return m_yvalues[0];
    if (coordinate>m_upper) return m_yvalues[m_nbin-1];
    for (int i=1;i<m_nbin;i++) {
      if ( (coordinate >= m_lower + (i-1)*m_binsize) &&
	   (coordinate <  m_lower + i*m_binsize) ) 
	return m_yvalues[i];
    }
  }
  return -1.0;
}

double Histogram::BinOrInterpolate(int bin) const {
  double binheight=m_yvalues[bin];
  if (binheight>0.0) {
    return binheight;
  }
  else {
    double leftbinheight(0.0), rightbinheight(0.0);
    for (int ibin=bin; ibin>0; --ibin) {
      leftbinheight=Bin(ibin-1);
      if (leftbinheight>0.0) break;
    }
    for (int ibin=bin; ibin<m_nbin-1; ++ibin) {
      rightbinheight=Bin(ibin+1);
      if (rightbinheight>0.0) break;
    }
    return (leftbinheight+rightbinheight)/2.0;
  }
}


double Histogram::BinOrInterpolate(double coordinate) const {
  if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

  int bin = int((coordinate-m_lower)/m_binsize+1.);
  if (bin<0) bin=0;
  if (bin>=m_nbin) bin=m_nbin-1;

  return BinOrInterpolate(bin);
}

void Histogram::Extrapolate(double coordinate,double * res,int mode) {
  if (!m_active) {
    msg_Error()<<"Error in Histogram : Tried to access a histogram with binsize <= 0 ! Return 0.."<<std::endl;
  }
  else {
    if (m_logarithmic>0) coordinate = log(coordinate)/m_logbase;

    for (int i=1;i<m_nbin;i++) {
      if ( (coordinate >= m_lower + (i-1)*m_binsize) &&
	   (coordinate <  m_lower + i*m_binsize) ) {
	res[0] = m_yvalues[i-1] +
	  (m_yvalues[i]-m_yvalues[i-1])/m_binsize *
	  (coordinate - m_lower - (i-1) * m_binsize);
	double over = 0.; double under = 0.;
	switch (mode) {
	case 1  : 
	  under = (coordinate - m_lower - (i-1) * m_binsize)/m_binsize * m_yvalues[i];
	  over = (m_lower + (i) * m_binsize - coordinate)/m_binsize * m_yvalues[i];
	  for (int j=0;j<i;j++) {
	    under += m_yvalues[j];
	  }
	  for (int j=i;j<m_nbin-1;j++) {
	    over += m_yvalues[j+1];
	  }
	  res[0]  = over;

	  if (m_depth>1) {
	    res[1]=0.;
	    for (int j=i;j<m_nbin;j++) {
	      res[1] = Max(res[1],m_y2values[j]);
	    }
	  }


	  break;
	case -1 :
	  for (int j=i-1;j>0;j--) {
	    over  += m_yvalues[j];
	    under += m_yvalues[j-1];
	  }
	  over    += m_yvalues[0];
	  res[0]  += (over+under)/2;
	  break;
	default : break;
	}
	/*
	if (m_depth>0) {
	  for (int j=1;j<=m_depth;j++) 
	    res[j] = Max(m_bins[i][j],m_bins[i-1][j]);
	}
	*/

      }
    }
  }
}

double Histogram::Integral() const 
{
  double total=0.;
  for (int i=0;i<m_nbin;i++) { 
    total += m_yvalues[i]; 
  }
  return total*m_binsize;
}

double Histogram::Ymax() const
{
  double ymax=m_yvalues[1];
  for (int i=1;i<m_nbin-1;i++) { 
    if (ymax<m_yvalues[i]) ymax=m_yvalues[i];
  }
  return ymax;
}

double Histogram::Ymin() const
{
  double ymin=1.e+65;
  for (int i=1;i<m_nbin-1;i++) { 
    if (ymin>m_yvalues[i] && m_yvalues[i]!=0) ymin=m_yvalues[i];
  }
  return ymin;
}

double Histogram::LogCoeff() const
{
  double ymax=m_yvalues[1];
  double ymin=1.e+65;
  double meany,meany2,meanly,meanly2;
  meany=meany2=meanly=meanly2=0.;

  int nl = 0;
  for (int i=1;i<m_nbin-1;i++) { 
    if (ymax<m_yvalues[i]) ymax=m_yvalues[i];
    if (ymin>m_yvalues[i] && m_yvalues[i]!=0.) ymin=m_yvalues[i];
    double y=m_yvalues[i];
    if (y!=0.) {
      meany   += y;
      meany2  += y*y;
      meanly  += log(y);
      meanly2 += sqr(log(y));
      ++nl;
    }
  }
  double rl = 0.;
  double rn = 0.;
  if (ymax!=0. && ymin!=0. && nl!=0) {
    double lymax=log(ymax);
    double lymin=log(ymin);
    meanly  = meanly/nl;
    meanly2 = meanly2/nl;
    double sigl2=meanly2-sqr(meanly);
    double ly0=0.5*(lymax+lymin);
    if (sigl2!=0.) rl=sigl2/sqr(ly0-meanly);
  }
  if (nl!=0) {
    meany   = meany/nl;
    meany2  = meany2/nl;
    double sig2=meany2-sqr(meany);
    double y0=0.5*(ymax+ymin);
    if (sig2!=0) rn=sig2/sqr(y0-meany);
  }
  double r;
  if (rl==0. && rn==0.) r=1.;
  else if (rl==0.)      r=0.;
  else if (rn==0.)      r=20.;
  else r=rl/rn;
  return r;
}



Histogram & Histogram::operator+=(const Histogram & histo)
{
  if (histo.m_nbin!=m_nbin) {
    msg_Error()<<"Error in Histogram : can not add histograms with different number of bins"<<std::endl;
    return *this;
  }
  for (int i=0;i<m_nbin;i++) { 
    m_yvalues[i]+= histo.m_yvalues[i]; 
  }
  if (m_depth>1) {
    for (int i=0;i<m_nbin;i++) { 
      m_y2values[i]+= histo.m_y2values[i]; 
    }
  }
  if (m_depth>2) {
    for (int i=0;i<m_nbin;i++) { 
      m_psvalues[i]+= histo.m_psvalues[i]; 
    }
  }
  
  m_fills+=histo.m_fills;
  m_psfills+=histo.m_psfills;
  return *this;
}


Histogram & Histogram::operator=(const Histogram & histo)
{
  for (int i=0;i<m_nbin;i++) { 
    m_yvalues[i]= histo.m_yvalues[i]; 
  }
  if (m_depth>1) {
    for (int i=0;i<m_nbin;i++) { 
      m_y2values[i]= histo.m_y2values[i]; 
    }
  }
  if (m_depth>2) {
    for (int i=0;i<m_nbin;i++) { 
      m_psvalues[i]= histo.m_psvalues[i]; 
    }
  }
  
  m_fills=histo.m_fills;
  m_psfills=histo.m_psfills;
  return *this;
}


void Histogram::Addopt(const Histogram & histo)
{
  
  if (m_depth<=1) {
    msg_Error()<<"Error in Histogram : can not Addopt histograms without statistical errors"<<std::endl;
    return;
  }
  if (histo.m_nbin!=m_nbin) {
    msg_Error()<<"Error in Histogram : can not add histograms with different number of bins"<<std::endl;
    return;
  }
  for (int i=0;i<m_nbin;i++) { 
    double w1=sqr(m_yvalues[i])/m_y2values[i];
    double w2=sqr(histo.m_yvalues[i])/histo.m_y2values[i];
    if (!(w1>0. && w2>0.)) w1=w2=1.;
    m_yvalues[i]=(m_yvalues[i]*w1+histo.m_yvalues[i]*w2)/(w1+w2); 
    m_y2values[i]=sqr(m_yvalues[i])/(w1+w2); 

    if (m_depth>2) {
      m_psvalues[i]+= histo.m_psvalues[i];       
    }
  }  
  m_fills+=histo.m_fills;
  m_psfills+=histo.m_psfills;
  return;
}

void Histogram::AddGeometric(const Histogram & histo)
{
  if (histo.m_nbin!=m_nbin) {
    msg_Error()<<"Error in Histogram : can not add histograms with different number of bins"<<std::endl;
    return;
  }
  for (int i=0;i<m_nbin;i++) {
    m_yvalues[i]=sqrt(BinOrInterpolate(i)*histo.BinOrInterpolate(i));
    if (m_depth>1 && histo.m_depth>1)
      m_y2values[i]=sqrt(m_y2values[i]*histo.m_y2values[i]);
  }
  m_fills+=histo.m_fills;
  m_psfills+=histo.m_psfills;
  return;
}

void Histogram::BinMin(const Histogram & histo)
{
  if (histo.m_nbin!=m_nbin) {
    msg_Error()<<"Error in Histogram::Min : histograms have different number of bins"<<std::endl;
    return;
  }
  for (int i=0;i<m_nbin;i++) { 
    double h1=m_yvalues[i];
    double h2=histo.m_yvalues[i];
    m_yvalues[i] = Min(h1,h2);
    if (m_depth>1) {
      if (h2<h1) m_y2values[i]= histo.m_y2values[i];
    }
    if (m_depth>2) {
      if (h2<h1) m_psvalues[i]= histo.m_psvalues[i];       
    }
  }  
  return;
}

void Histogram::BinMax(const Histogram & histo)
{
  if (histo.m_nbin!=m_nbin) {
    msg_Error()<<"Error in Histogram::Max : histograms have different number of bins"<<std::endl;
    return;
  }
  for (int i=0;i<m_nbin;i++) { 
    double h1=m_yvalues[i];
    double h2=histo.m_yvalues[i];
    m_yvalues[i] = Max(h1,h2);
    if (m_depth>1) {
      if (h2>h1) m_y2values[i]= histo.m_y2values[i];
    }
    if (m_depth>2) {
      if (h2>h1) m_psvalues[i]= histo.m_psvalues[i];       
    }
  }  
  return;
}

int Histogram::CheckStatistics(const Histogram & histo,double& avgs,double& stdd)
{
  if (!m_finished||(!(histo.m_finished))) {
    msg_Error()<<"Error in Histogram : Histogram must be Finalized for CheckStatistics()!"<<std::endl;
    return 0;
  }
  if (m_depth<=1) {
    msg_Error()<<"Error in Histogram : can not CheckStatistics() histograms without statistical errors"<<std::endl;
    return 0;
  }
  avgs=0.;stdd=0.;
  double cnt=0.;
  for (int i=1;i<m_nbin-1;i++) { 
    double dxs=0.;
    if (!IsEqual(m_y2values[i],sqr(m_yvalues[i]))&&
	!IsEqual(histo.m_y2values[i],sqr(histo.m_yvalues[i]))) {
      dxs=(m_yvalues[i]-histo.m_yvalues[i])/sqrt(m_y2values[i]+histo.m_y2values[i]);
      cnt+=1.;
    }
    avgs+=dxs;
    stdd+=sqr(dxs);
  }
  if (cnt>0.) {
    avgs/=cnt;
    stdd/=cnt;
    stdd=sqrt(stdd);
  }
  return int(cnt);
}