#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 "FOpenCalibFile.c"

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


vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;

unsigned int CurrentEventID;

bool breakout=false;

size_t ROIxNOP;
UInt_t NumberOfPixels;
UInt_t RegionOfInterest;
int NEvents;

size_t drs_n;
vector<float> drs_basemean;
vector<float> drs_gainmean;
vector<float> drs_triggeroffsetmean;

int FOpenDataFile( fits & );


vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
vector<float> Vdiff(FAD_MAX_SAMPLES);  // numerical derivative

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


float getValue( int, int );
void computeN1mean( int );
void removeSpikes( int );

// 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

TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
TH1F *hMeanBsl, *hpltMeanBsl;
TH1F *hRmsBsl, *hpltRmsBsl; 
TH2F * hAmplSpek_cfd;
TH2F * hAmplSpek_discri;
TH2F * hTemp_Array[1441];
TObjArray hList;
TObjArray hListBaseline, hListTemplates;

void BookHistos( );
void SaveHistograms( char * );

int searchSinglesPeaks = 0;

int fpeak(
	char *datafilename 		= "../../20111011_055.fits.gz",
	const char *drsfilename = "../../20111011_054.drs.fits.gz",
	int nevents 			= -1,
	int firstevent 			= 0,
	int PixelID				= -1,
	bool spikeDebug = false,
	int verbosityLevel = 1 // different verbosity levels can be implemented here
	)
{

// Create (pointer to) Canvases, which are used in every run,
// also in 'non-debug' runs
	TCanvas * cSpektrum;
	TCanvas * cStartCell;
	TCanvas *cTemplate;
	cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
	cSpektrum->Divide(1,2);
	cStartCell = new TCanvas("cStartCell ", "The Startcells of this run", 10,410,400,400);
	cTemplate = new TCanvas("cTemplate","Template of current Pixel",1,1,1600,1000);

	// Canvases only need if spike Debug, but I want to deklare
	// the pointers anyway ...
	TCanvas *cRawAndSpikeRemoval = NULL;
	TCanvas *cFiltered = NULL;


	if (spikeDebug){
		cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",410,10,400,400);
		cRawAndSpikeRemoval->Divide(1, 2);
		cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
		cFiltered->Divide(1, 2);
	}

gStyle->SetPalette(1,0);
gROOT->SetStyle("Plain");
// read FACT raw data 
// 	* remove spikes
//	* calculate baseline
//	* find peaks (CFD and normal discriminator)
//	* compute pulse height and pulse integral spektrum of the peaks

// sliding window filter settings
//	int k_slide = 16;
//	vector<double> a_slide(k_slide, 1);
//	double b_slide = k_slide;

// 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.;

// 2nd slinding window filter
//	int ks2 = 16;
//	vector<double> as2(ks2, 1);
//	double bs2 = ks2;
	
// Open the data file
	fits *datafile = new fits( datafilename );
	if (!datafile) {
	    printf( "Could not open the file: %s\n", datafilename );
		  return 1;
	}
	
// access data
	NEvents = FOpenDataFile( *datafile );
	printf("number of events in file: %d\n", NEvents);
  if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!

//Get the DRS calibration
	FOpenCalibFile(	drsfilename,
					drs_basemean,
					drs_gainmean,
					drs_triggeroffsetmean,
					drs_n);

//Check the sizes of the data columns
	if (drs_n != 1474560) {
		cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
			<< drs_n << endl;
		cerr << " Aborting." << endl;
		return 1;
	}

	if(ROIxNOP != 1474560)
	{
		cout << "warning: data_n should better be 1440x1024=1474560, but is "
			<< ROIxNOP << endl;
		cout << "this script is not guaranteed to run under these "
			<<" circumstances....any way ... it is never guaranteed." << endl;
	}
// Book the histograms

	BookHistos( );

	float calibratedVoltage;

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

		for ( int pix = 0; pix < 1440; pix++ ){

			hStartCell->Fill( pix, StartCellVector[pix] );

			// this is a stupid hack ... there is more code at the
			// end of this loop to complete this hack ...
			// beginning with if (PixelID != -1) as well.
			if (PixelID != -1) {
				pix = PixelID;
				if (verbosityLevel > 0){
					cout << "Processing Event number: " << CurrentEventID << "\t"
						<< "Pixel number: "<< pix << endl;
				}
			}

			if (verbosityLevel > 0){
				if (pix % 20 ==0){
					cout << "Processing Event number: " << CurrentEventID << "\t"
						<< "Pixel number: "<< pix << endl;
				}
			}

			// compute the DRs calibrated values and put them into Ameas[]
			for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
				calibratedVoltage = getValue( sl, pix);
				if (verbosityLevel > 10){
					printf("calibratedVoltage = %f\n", calibratedVoltage);
				}
				Ameas[ sl ] = calibratedVoltage;

				// in case we want to plot this ... we need to put it into a
				// Histgramm
				if (spikeDebug){
					debugHistos[Ameas_].SetBinContent(sl, calibratedVoltage);
				}
			}
			// operates on Ameas[] and writes to N1mean[]
			computeN1mean( RegionOfInterest );

			// operates on Ameas[] and N1mean[], then writes to Vcorr[]
			removeSpikes( RegionOfInterest );

			// filter Vcorr with sliding average using FIR filter function
			//factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
			sliding_avg(Vcorr, Vslide, 8);
			// 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
			//factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
			sliding_avg(Vcfd, Vcfd2, 8);


			// 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"
			//ShiftRegionBy(*zXings, -ks2/2);
			EnlargeRegion(*zXings, 10, 10);
			findAbsMaxInRegions(*zXings, Vslide);
			removeMaximaBelow( *zXings, 11.0, 0);
			removeMaximaAbove( *zXings, 14.0, 0);

			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);

					for (int j=-100; j<150; j++){
						if (zXings->at(i).maxPos + j >= 0 && zXings->at(i).maxPos + j <= 1023)
							hTemp_Array[pix]->Fill(j+101, Vslide[zXings->at(i).maxPos + j]);
					}
				}
			}

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

				cRawAndSpikeRemoval->cd( 1);
				gPad->SetGrid();
				debugHistos[Ameas_].Draw();

				cRawAndSpikeRemoval->cd( 2);
				gPad->SetGrid();
				debugHistos[Vcorr_].Draw();

				cFiltered->cd(1);
				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(2);
				gPad->SetGrid();
				debugHistos[Vcfd2_].Draw();
				TLine *zeroline = new TLine(0, 0, 1024, 0);
				zeroline->SetLineColor(kBlue);
				zeroline->Draw();

				cRawAndSpikeRemoval->Update();
				cFiltered->Update();

				cTemplate->cd();
				hTemp_Array[pix]->Draw("COLZ");
				cTemplate->Modified();
				cTemplate->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 ...
				// all the Objects on the heap should be deleted here.

			}// end of if(spikeDebug)

			delete zXings;

			// this is the 2nd part of the ugly hack in the beginning.
			// this makes sure, that the loop ends
			if (PixelID != -1){
				pix = 1440;
			}


			if (breakout)
				break;

		} // end of loop over pixels

if (ev % 10 == 0){
		cSpektrum->cd(1);
		hAmplSpek_cfd->Draw("COLZ");
		cSpektrum->cd(2);
		hAmplSpek_discri->Draw("COLZ");
		cSpektrum->Modified();
		cSpektrum->Update();

		// updating seems not to work ..
		// debug cout
		cStartCell->cd();
		hStartCell->Draw();
		cStartCell->Modified();
		cStartCell->Update();

		cTemplate->cd();
		hTemp_Array[PixelID]->GetYaxis()->SetRangeUser(-10,40);
		hTemp_Array[PixelID]->Draw("COLZ");
		cTemplate->Modified();
		cTemplate->Update();
}


		if (breakout)
			break;
	}	// end of loop over pixels


	SaveHistograms( datafilename );

	return( 0 );
}

void removeSpikes(int Samples){

    const float fract = 0.8;
    float x, xp, xpp, x3p;

    // assume that there are no spikes
    for ( int i = 0; i <  Samples; i++) Vcorr[i] = Ameas[i];

// find the spike and replace it by mean value of neighbors
    for ( int i = 0; i < Samples; i++) {

    // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );

    x = Ameas[i] - N1mean[i];

        if ( x < -5. ){ // a spike candidate
            // check consistency with a single channel spike
            xp = Ameas[i+1] - N1mean[i+1];
            xpp = Ameas[i+2] - N1mean[i+2];
            x3p = Ameas[i+3] - N1mean[i+3];

            // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);

            if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
                // printf("double spike candidate\n");
                Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
                Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
                // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
                // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
                i = i + 3; 
            }
            else{
    
                if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
                    Vcorr[i+1] = N1mean[i+1];
                    // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
                    // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
                    N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
                    i = i + 2;//do not care about the next sample it was the spike
                }
                // treatment for the end of the pipeline must be added !!!
            }
        }
        else{
             // do nothing
        }
    } // end of spike search and correction
	for ( int i = 0; i < Samples; i++ ) debugHistos[ Vcorr_ ].SetBinContent( i, Vcorr[i] );
}

void computeN1mean( int Samples ){
// compute the mean of the left and right neighbors of a channel

    for( int i = 2; i < Samples - 2; i++){
/*        if (i == 0){ // use right sample es mean
            N1mean[i] = Ameas[i+1];
        }
        else if ( i == Samples-1 ){ //use left sample as mean
            N1mean[i] = Ameas[i-1];
        }
        else{
            N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
        }
*/
        N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
    }
} // end of computeN1mean computation

float getValue( int slice, int pixel ){
	const float dconv = 2000/4096.0;

	float vraw, vcal;

	unsigned int pixel_pt;
	unsigned int slice_pt;
	unsigned int cal_pt;
	unsigned int drs_cal_offset;

	// printf("pixel = %d, slice = %d\n", slice, pixel);

	pixel_pt = pixel * RegionOfInterest;
	slice_pt = pixel_pt + slice;
	drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
	cal_pt    = pixel_pt + drs_cal_offset;

	vraw = AllPixelDataVector[ slice_pt ] * dconv;
	vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;

	return( vcal );
}

// booking and parameter settings for all histos
void BookHistos( ){
	// histograms for baseline extraction
	char hName[500];
	char hTitle[500];

	TH1F *h;

	for( int i = 0; i < NPIX; i++ ) {
		sprintf(&hTitle[0],"all events all slices of pixel %d", i);
		sprintf(&hName[0],"base%d", i);

		h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );

		h->GetXaxis()->SetTitle( "Sample value (mV)" );
		h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
		hListBaseline.Add( h );
		hBaseline[i] = h;
	}

    hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
    hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
    hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    hList.Add( hMeanBsl );

    hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
    hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
    hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
    hList.Add( hpltMeanBsl );

    hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
    hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
    hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    hList.Add( hRmsBsl );

    hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
    hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
    hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
    hList.Add( hpltRmsBsl );

	hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",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 );

	hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
	hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
	hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
	hList.Add( hAmplSpek_discri );

	hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
	hStartCell->GetXaxis()->SetTitle( "pixel" );
	hStartCell->GetXaxis()->SetTitle( "slice" );
	hList.Add( hStartCell );

	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.)");
	}

	TH2F *temp;
	for ( int type = 0; type < 1441; type++){
		sprintf(&hTitle[0],"pulse template of pixel %d", type );
		sprintf(&hName[0],"template_%d", type );

		temp = new TH2F( hName, hTitle, 256, 0, 255, 256, -10.5, 117.5);

		temp->GetXaxis()->SetTitle( "Time slice (%.1f ns/slice)" );
		temp->GetYaxis()->SetTitle( "Amplitude (about mV)" );
		hListTemplates.Add( temp );
		hTemp_Array[ type ] = temp;
	}


}

void SaveHistograms( char * loc_fname ){

	TString fName; // name of the histogram file

	// create the filename for the histogram file
	fName = loc_fname; // use the name of the tree file
	// TODO ... next statement doesn't work for ".fits.gz"
	fName.Remove(fName.Length() - 5); // remove the extension .fits
	fName += "_discri.root"; // add the new extension

	// create the histogram file (replace if already existing)
	TFile tf( fName, "RECREATE");

	hList.Write(); // write the major histograms into the top level directory
	tf.mkdir("BaselineHisto");
	tf.cd("BaselineHisto"); // go to new subdirectory
	hListBaseline.Write(); // write histos into subdirectory
	tf.cd("..");
	tf.mkdir("PulseTempplates");
	tf.cd("PulseTempplates");
	hListTemplates.Write(); // write histos into subdirectory

	tf.Close(); // close the file
} // end of SaveHistograms( char * loc_fname )

int FOpenDataFile(fits &datafile){

//	cout << "-------------------- Data Header -------------------" << endl;
//	datafile.PrintKeys();
//	cout << "------------------- Data Columns -------------------" << endl;
//	datafile.PrintColumns();

	//Get the size of the data column
	RegionOfInterest	= datafile.GetUInt("NROI");
	NumberOfPixels		= datafile.GetUInt("NPIX");
	//Size of column "Data" = #Pixel x ROI
	ROIxNOP				= datafile.GetN("Data");

	//Set the sizes of the data vectors
	AllPixelDataVector.resize(ROIxNOP,0);
	StartCellVector.resize(NumberOfPixels,0);

	//Link the data to variables
	datafile.SetRefAddress("EventNum",  CurrentEventID);
	datafile.SetVecAddress("Data", AllPixelDataVector);
	datafile.SetVecAddress("StartCellData", StartCellVector);
	datafile.GetRow(0);

	return datafile.GetNumRows() ;
}