#include <iostream>
#include <fstream>
#include <string>
#include <vector>

#include <math.h> // exp

#define TOT_SAMPLES 8192
#define B_SAMPLES 400


#define GAIN 6.22575 //1/(POW(POLES,POLES)*EXP(-POLES)/POLES!)
#define T_C 100.0 			// Time constant (RC=CR)
#define PZ 500.0 			// Pole zero correction

void Filter(std::vector<double>* wave) // Semi-gaussian filter
{	
	
	//std::vector<double> old_wave = *wave;
	double b1,a0,a1,a0_1;
	b1=exp(-1./T_C);
	a0=(1.+b1)/2;
	a0_1=1.-b1;

	int i,j;
	
	// 1 differentiator (high-pass filter)
	double my_value, baseline;
	
	baseline= (*wave)[0];
	(*wave)[0]=0;
	for(i=1;i<wave->size();++i)
	{
		my_value=(*wave)[i]-baseline; // Baseline
		(*wave)[i]=b1*(*wave)[i-1] + a0*(my_value-(*wave)[i-1]) + (*wave)[i-1]/PZ;
	}
	
	//  6 integrators (low-pass filter)
 	for(j=0;j<6;++j)
	{
		for(i=1;i<wave->size();++i) 
		{
			my_value=b1*(*wave)[i-1]+a0_1*(*wave)[i];
			(*wave)[i]=my_value;
		}
	}
	
	//for(i=0;i<wave->size();++i) std::cout<<old_wave[i]-old_wave[0]<<" "<<(*wave)[i]<<std::endl;
};

unsigned int GetAmplitude(std::vector<double> wave)
{
	double minimum=wave[0];
	double baseline=0;
	int i;	
	
	for(i=0;i<B_SAMPLES && i<wave.size();++i) baseline+=wave[i];
	baseline*=1./B_SAMPLES;

	for(;i<wave.size();++i) if(wave[i]<minimum) minimum=wave[i];

	return baseline-minimum;

};

int main (int argc, char* argv[])
{
	std::fstream f_in;
	std::string str;
	if(argc<2)
	{
		std::cerr<<"USAGE: ./WavesToHisto input_file"<<std::endl;
		return 1;
	}
	str=argv[1];
	f_in.open(str.c_str(),std::fstream::in);
	if(!f_in)
	{
		std::cerr<<"No such file: "<<str<<std::endl;
	}

	// Header
	std::getline(f_in,str);
	while(str!="")
	{
		std::getline(f_in,str);
	}

	// Waveforms
	int read;
	int i=0;
	double baseline, minimum;

	int histo[500]={0};
	unsigned int channel;
	std::vector<double> wave;

	while(f_in)
	{
		i=0;
		baseline=0;
		f_in>>read;
		while( f_in && i<(TOT_SAMPLES-3) )
		{
			wave.push_back(read);
			f_in>>read;
			++i;
		}
		f_in>>read>>read;
		if(wave.size()>0)
		{		
			Filter(&wave);
			channel=GetAmplitude(wave);
			if(channel<500) histo[channel]++;
		}

		wave.clear();
				

		/*		
		minimum=read;	
		while(i<(TOT_SAMPLES-3))
		{
			//std::cout<<i<<" "<<read<<std::endl;
			if(i<B_SAMPLES) baseline+=read;
			else if (read<minimum) minimum=read;

			f_in>>read;
			i++;
		}
		f_in>>read>>read; // Waveform tail
		baseline*=1./B_SAMPLES;
		channel=baseline-minimum;
		if(channel<500) histo[channel]++;
		//std::cout<<baseline<<" "<<minimum<<" "<<baseline-minimum<<std::endl;
		*/
	}

	for(i=0;i<500;++i) std::cout<<i<<" "<<histo[i]<<std::endl;
	
	return 0;
}