// $Id$
//
//    File: DHistogram.h
// Created: Mon Aug  1 08:56:58 EDT 2011
// Creator: davidl (on Linux ifarm1102 2.6.18-128.7.1.el5 x86_64)
//

#ifndef _DHistogram_
#define _DHistogram_


#include <iostream>
#include <cstddef> // for NULL
#include <cmath>

#include <TH1.h>

#include <JANA/jerror.h>

class DHistogram{
	public:
		inline DHistogram(int Nbins, double lowEdge, double highEdge);
		inline virtual ~DHistogram();
		
		inline void Add(const DHistogram *h);
		
		inline double Fill(double x, double weight=1.0);
		inline int FindBin(double x) const;
		inline int FindFirstBinAbove(double val) const;
		inline int FindFirstNonZeroBin(void) const;
		inline int FindLastNonZeroBin(void) const;
		
		inline int GetNbins(void) const;
		inline double GetBinLowEdge(int bin) const;
		inline double GetBinCenter(int bin) const;
		inline double GetBinContent(int bin) const;
		inline double GetBinWidth(void) const;
		
		inline double Integral(void) const;
		inline TH1D* MakeTH1D(string name, string title) const;
		
		inline void Reset(void);
		inline void Scale(double scale);
		inline void SetBinContent(int bin, double val);
		
	protected:
	
	
	private:
		DHistogram(); // prevent default constructor usage
	
		int Nbins;
		double lowEdge;
		double highEdge;
		double binWidth;
		double half_binWidth;

		double *content;
		double *xvals; // bin center
};


//---------------------------------
// DHistogram    (Constructor)
//---------------------------------
inline DHistogram::DHistogram(int Nbins, double lowEdge, double highEdge)
{
	
	this->Nbins = Nbins;
	this->lowEdge = lowEdge;
	this->highEdge = highEdge;
	binWidth = (highEdge - lowEdge)/(double)Nbins;
	half_binWidth = binWidth/2.0;

	if(Nbins>0 && Nbins<100000000){
		content = new double[Nbins+1];
		xvals = new double[Nbins+1];
		
		for(int bin=1; bin<=Nbins; bin++){
			content[bin] = 0.0;
			xvals[bin] = lowEdge + (double)(bin-1)*binWidth + half_binWidth;
		}
	}else{
		content = NULL;
		xvals = NULL;
	}
}

//---------------------------------
// ~DHistogram    (Destructor)
//---------------------------------
inline DHistogram::~DHistogram()
{
	if(content!=NULL) delete[] content;
	if(xvals!=NULL) delete[] xvals;
}

//---------------------------------
// Add
//---------------------------------
inline void DHistogram::Add(const DHistogram *h)
{
	if(h->GetBinWidth()!=binWidth || h->GetNbins()!=Nbins){
		_DBG_<<"Histogram parameters don't match!!"<<endl;
		exit(-1);
	}
	
	for(int bin=1; bin<=Nbins; bin++){
		content[bin] += h->GetBinContent(bin);
	}
}

//---------------------------------
// FindBin
//---------------------------------
inline int DHistogram::FindBin(double x) const
{
	int bin = 1 + (int)floor((x-lowEdge)/binWidth);
	
	return bin;
}

//---------------------------------
// FindFirstBinAbove
//---------------------------------
inline int DHistogram::FindFirstBinAbove(double val) const
{
	for(int bin=1; bin<=Nbins; bin++){
		if(content[bin] >= val)return bin;
	}
	
	return -1;
}

//---------------------------------
// FindFirstNonZeroBin
//---------------------------------
int DHistogram::FindFirstNonZeroBin(void) const
{
	for(int bin=1; bin<=Nbins; bin++){
		if(content[bin] != 0.0)return bin;
	}
	
	return -1;
}

//---------------------------------
// FindLastNonZeroBin
//---------------------------------
int DHistogram::FindLastNonZeroBin(void) const
{
	for(int bin=Nbins; bin>=1; bin--){
		if(content[bin] != 0.0)return bin;
	}
	
	return -1;
}

//---------------------------------
// Fill
//---------------------------------
inline double DHistogram::Fill(double x, double weight)
{
	int bin = 1 + (int)floor((x-lowEdge)/binWidth);
	if(bin<1 || bin>Nbins)return 0.0;

	content[bin] += weight;
	
	return content[bin];
}
		
//---------------------------------
// GetNbins
//---------------------------------
inline int DHistogram::GetNbins(void) const
{
	return Nbins;
}

//---------------------------------
// GetBinLowEdge
//---------------------------------
inline double DHistogram::GetBinLowEdge(int bin) const
{
	return lowEdge + (double)(bin-1)*binWidth;
}

//---------------------------------
// GetBinCenter
//---------------------------------
inline double DHistogram::GetBinCenter(int bin) const
{
	if(bin<1 || bin>Nbins)return 0.0;

	return xvals[bin];
}

//---------------------------------
// GetBinContent
//---------------------------------
inline double DHistogram::GetBinContent(int bin) const
{
	if(bin<1 || bin>Nbins)return 0.0;

	return content[bin];
}

//---------------------------------
// GetBinWidth
//---------------------------------
inline double DHistogram::GetBinWidth(void) const
{
	return binWidth;
}

//---------------------------------
// Integral
//---------------------------------
inline double DHistogram::Integral(void) const
{
	double I=0.0;

	for(int bin=1; bin<=Nbins; bin++) I += content[bin];

	return I;
}

//---------------------------------
// MakeTH1D
//---------------------------------
inline TH1D* DHistogram::MakeTH1D(string name, string title) const
{
	TH1D *h = new TH1D(name.c_str(), title.c_str(), Nbins, lowEdge, highEdge);

	for(int bin=1; bin<=Nbins; bin++) h->SetBinContent(bin, content[bin]);

	return h;
}

//---------------------------------
// Reset
//---------------------------------
inline void DHistogram::Reset(void)
{
	for(int bin=1; bin<=Nbins; bin++)content[bin] = 0.0;
}

//---------------------------------
// Scale
//---------------------------------
inline void DHistogram::Scale(double scale)
{
	for(int bin=1; bin<=Nbins; bin++)content[bin] *= scale;
}

//---------------------------------
// SetBinContent
//---------------------------------
inline void DHistogram::SetBinContent(int bin, double val)
{
	if(bin<1 || bin>Nbins)return;
	
	content[bin] = val;
}


#endif // _DHistogram_