#include <TROOT.h>
#include <TCanvas.h>
#include <TProfile.h>
#include <TTimer.h>
#include <TH1F.h>
#include <TH2F.h>
#include <Getline.h>
#include <TLine.h>
#include <TBox.h>
#include <TMath.h>
#include <TFile.h>
#include <TStyle.h>

#include <stdio.h>
#include <stdint.h>
#include <cstdio>
#include <deque>

#define NPIX  1440
#define NCELL 1024
#define FAD_MAX_SAMPLES 1024

#define HAVE_ZLIB
#include "fits.h"

#include "openFits.h"
#include "openFits.c"

#include "discriminator.h"
#include "discriminator.C"
#include "zerosearch.h"
#include "zerosearch.C"
#include "factfir.C"

#include "FOpenDataFile.h"
#include "FOpenDataFile.c"

#include "DrsCalibration.C"
#include "DrsCalibration.h"

#include "SpikeRemoval.h"
#include "SpikeRemoval.C"

bool breakout=false;

int NEvents;
vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;
unsigned int CurrentEventID;
size_t PXLxROI;
UInt_t RegionOfInterest;
UInt_t NumberOfPixels;

size_t TriggerOffsetROI, RC;
vector<float> Offset, Gain, TriggerOffset;

vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result
vector<float> Vcfd(FAD_MAX_SAMPLES);    // CDF result
vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average

// histograms
const int NumberOfDebugHistoTypes = 7;
const unsigned int
	Ameas_ 	= 0,
	N1mean_ = 1,
	Vcorr_  = 2,
	Vtest_  = 3,
	Vslide_ = 4,
	Vcfd_   = 5,
	Vcfd2_  = 6;

TH1F* h;
TH1F *debugHistos;
TH2F* hStartCell;		// id of the DRS physical pipeline cell where readout starts
						// x = pixel id, y = DRS cell id

TH2F * hAmplSpek_cfd;
TObjArray hList;

void BookHistos( TString &title );
void SaveHistograms(const char * );

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// read FACT raw data
// 	* remove spikes
//	* calculate baseline
//	* find peaks (CFD and normal discriminator)
//	* compute pulse height and pulse integral spektrum of the peaks
int fpeak_cfd(
	char *datafilename 		= "data/20111016_013.fits.gz",
	const char *drsfilename = "../../20111016_011.drs.fits.gz",
	const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root",
	int firstevent 			= 0,
	int nevents 			= -1,
	int firstpixel			= 0,
	int npixel				= -1,
	bool spikeDebug = false,
	int avg1		= 8,
	int avg2 		= 8,
	float threshold = 5.0,
	int verbosityLevel = 1, // different verbosity levels can be implemented here
	bool ProduceGraphic = false
	)
{
gStyle->SetPalette(1,0);
gROOT->SetStyle("Plain");

TString histogramtitle;
histogramtitle += datafilename; 
histogramtitle += " ";
histogramtitle += drsfilename ;
histogramtitle += " ";

histogramtitle += avg1;
histogramtitle += " ";
histogramtitle +=  avg2;
histogramtitle += " ";
histogramtitle +=  threshold ;
histogramtitle +=  " 100";
cout << histogramtitle.Data() << endl;


	TCanvas *cFiltered = NULL;
	if (spikeDebug){
		cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
		cFiltered->Divide(1, 4);
	}



// CFD filter settings
	int k_cfd = 10;
	vector<double> a_cfd(k_cfd, 0);
	double b_cfd = 1.;
	a_cfd[0]=-0.75;
	a_cfd[k_cfd-1]=1.;

	fits * datafile;
	// Opens the raw data file and 'binds' the variables given as
	// Parameters to the data file. So they are filled with
	// raw data as soon as datafile->GetRow(int) is called.
	NEvents = openDataFits( datafilename, &datafile,
		AllPixelDataVector, StartCellVector, CurrentEventID,
		RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
	if (NEvents == 0){
		cout << "return code of OpenDataFile:" << datafilename<< endl;
		cout << "is zero -> aborting." << endl;
		return 1;
	}

	if (verbosityLevel > 0)
		cout <<"number of events in file: "<< NEvents << "\t";
	if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
	if (verbosityLevel > 0)
		cout <<"of, which "<< nevents << "will be processed"<< endl;

	if (verbosityLevel > 0)
		cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
	if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
	if (verbosityLevel > 0)
		cout <<"of, which "<< npixel << "will be processed"<< endl;

	RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI);
	if (RC == 0){
		cout << "return code of openCalibFits:" << drsfilename << endl;
		cout << "is zero -> aborting." << endl;
		return 1;
	}

	BookHistos( histogramtitle  );
//-----------------------------------------------------------------------------
// Loops over Every Event and Pixel
//-----------------------------------------------------------------------------
	cout << "Processing Events in total: " << nevents << endl;

	for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
		// Get an Event --> consists of 1440 Pixel ...erm....data
		datafile->GetRow( ev );

		for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
			if (verbosityLevel > 0){
				if ( pix == (firstpixel + npixel - 1) ) {
					if (ev == firstevent) {
					    printf ("%06d done", CurrentEventID);
					    fflush(stdout);
					} else {
					    printf ("\b\b\b\b\b\b\b\b\b\b\b%06d done", CurrentEventID);
					    fflush(stdout);
					}
				}
			}

			if ((pix+1)%9==0){
			applyDrsCalibration( Ameas,pix,12,60,
				Offset, Gain, TriggerOffset,
				RegionOfInterest, AllPixelDataVector, StartCellVector);
			}else{
			applyDrsCalibration( Ameas,pix,12,12,
				Offset, Gain, TriggerOffset,
				RegionOfInterest, AllPixelDataVector, StartCellVector);
			}

			// finds spikes in the raw data, and interpolates the value
			// spikes are: 1 or 2 slice wide, positive non physical artifacts
			removeSpikes (Ameas, Vcorr);

			// filter Vcorr with sliding average using FIR filter function
			sliding_avg(Vcorr, Vslide, avg1);
			// filter Vslide with CFD using FIR filter function
			factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
			// filter Vcfd with sliding average using FIR filter function
			sliding_avg(Vcfd, Vcfd2, avg2);


			// peaks in Ameas[] are found by searching for zero crossings
			// in Vcfd2
			// first Argument 1 means ... *rising* edge
			// second Argument 1 means ... search with stepsize 1 ... 8 is okay as well
			vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
			// zXings means "zero cross ings"
			EnlargeRegion( *zXings, 10, 10);
			findAbsMaxInRegions( *zXings, Vslide);
			removeMaximaBelow( *zXings, threshold, 0);
			removeRegionWithMaxOnEdge( *zXings, 2);
			removeRegionOnFallingEdge( *zXings, 100);
			
			// fill maxima in Histogram
			if (zXings->size() != 0 ){
				for (unsigned int i=0; i<zXings->size(); i++){
					if (verbosityLevel > 1){
						cout << zXings->at(i).maxPos << ":\t"
							<< zXings->at(i).maxVal <<endl;
					}
					hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
				}
			}

			if ( spikeDebug ){
		
					// TODO do this correct. The vectors should be the rigt ones... this is just luck	
					debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
					debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
					debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
					debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
					debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);

				for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
					debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
					debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
					debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
					debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
					debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
				}

				cFiltered->cd(1);
				gPad->SetGrid();
				debugHistos[Ameas_].Draw();
				
				cFiltered->cd(2);
				gPad->SetGrid();
				debugHistos[Vcorr_].Draw();

				cFiltered->cd(3);
				gPad->SetGrid();
				debugHistos[Vslide_].Draw();

				TBox *OneBox;
				vector<TBox*> MyBoxes;
				for (unsigned int i=0; i<zXings->size(); i++){
					OneBox = new TBox(
						zXings->at(i).maxPos -10 ,
						zXings->at(i).maxVal -0.5,
						zXings->at(i).maxPos +10 ,
						zXings->at(i).maxVal +0.5);
					OneBox->SetLineColor(kBlue);
					OneBox->SetLineWidth(1);
					OneBox->SetFillStyle(0);
					OneBox->SetFillColor(kRed);
					MyBoxes.push_back(OneBox);
					OneBox->Draw();
				}

				cFiltered->cd(4);
				gPad->SetGrid();
				debugHistos[Vcfd2_].Draw();
				TLine *zeroline = new TLine(0, 0, 1024, 0);
				zeroline->SetLineColor(kBlue);
				zeroline->Draw();

				cFiltered->Update();
				//Process gui events asynchronously during input
				TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
				timer.TurnOn();
				TString input = Getline("Type 'q' to exit, <return> to go on: ");
				timer.TurnOff();
				if (input=="q\n") {
					breakout=true;
				}

				//TODO!!!!!!!!!
				// do some Garbage collection ...
			}// end of if(spikeDebug)

			delete zXings;
			if (breakout)
				break;

		} // end of loop over pixels

		if (breakout)
			break;
	}	// end of loop over pixels
	
	cout << "Event processing done." << endl;  

if (ProduceGraphic){
	TCanvas * cSpektrum;
	cSpektrum = new TCanvas("cSpektrum","Amplitude Spektrum",10,10,400,400);
	cSpektrum->Divide(1,1);
  cSpektrum->cd(1);
  hAmplSpek_cfd->Draw("COLZ");
  cSpektrum->Modified();
  cSpektrum->Update();
}
	if ( strlen(OutRootFileName) != 0){
		SaveHistograms( OutRootFileName );
	}
	return( 0 );
}

// booking and parameter settings for all histos
void BookHistos( TString& title){
	// histograms for baseline extraction

	

	hAmplSpek_cfd = new TH2F("hAmplSpek_cfd", title,1440,-0.5,1439.5, 256, -27.5, 100.5);
	hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
	hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
	hList.Add( hAmplSpek_cfd );


	debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
	for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
		debugHistos[ type ].SetBins(1024, 0, 1024);
		debugHistos[ type ].SetLineColor(1);
		debugHistos[ type ].SetLineWidth(2);

		// set X axis paras
		debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
		debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
		debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
		debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));

		// set Y axis paras
		debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
		debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
		debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
		debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
	}
}

void SaveHistograms( const char * loc_fname){
	// create the histogram file (replace if already existing)
	TFile tf( loc_fname, "RECREATE");

	hList.Write(); // write the major histograms into the top level directory

	tf.Close(); // close the file
}