#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 <TMath.h>

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


#define HAVE_ZLIB
#include "izstream.h"
#include "fits.h"
//#include "TPKplotevent.c"
//#include "FOpenDataFile.c"
#include "FOpenCalibFile.c"

#include "discriminator.h"
#include "discriminator.C"

#include "zerosearch.C"

#define NPIX  1440
#define NCELL 1024

// data access and treatment
#define FAD_MAX_SAMPLES 1024

vector<int16_t> data;
vector<int16_t> data_offset;

unsigned int data_num;

size_t data_n;
UInt_t data_px;
UInt_t data_roi;
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> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 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

#include "factfir.C"

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

// histograms
const int Ntypes = 7;
const unsigned int  // arranged by Dominik
	tAmeas 	= 0, 
	tN1mean = 1, 
	tVcorr  = 2,
	tVtest  = 3,
	tVslide = 4,
	tVcfd   = 5,
	tVcfd2  = 6;

TH1F* h;
TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
                    // x = pixel id, y = DRS cell id
TH2F hPixelCellData("PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);

void BookHistos( int );

// Create a canvas
TCanvas* CW;
TCanvas* cFilter;

int spikeDebug = 0;


int zippedfitstest( 
	char *datafilename 		= "../raw/20111011_026.fits.gz",
	const char *drsfilename = "../raw/20111011_022.drs.fits.gz",
	int pixelnr 			= 4,
	int firstevent 			= 0, 
	int nevents 			= -1 ){
// read and analyze FACT raw data

// 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;
	gROOT->SetStyle("Plain");
	
//-------------------------------------------
// Open the file
//-------------------------------------------
	fits datafile( 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);

    // compare the number of events in the data file with the nevents the
    // the user would like to read. nevents = -1 means: read all
    if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;

	
//-------------------------------------------
//Get the DRS calibration
//-------------------------------------------

	FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
	
//-------------------------------------------
//Check the sizes of the data columns
//-------------------------------------------
	if(drs_n!=data_n)
	{
		cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
		return 1;
	}
// Book the histograms
    BookHistos( data_roi );

// Loop over events
	cout << "--------------------- Data --------------------" << endl;

    float value;
	TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5);

    for ( int ev = firstevent; ev < firstevent + nevents; ev++) {

	    datafile.GetRow( ev );
		
		if (ev % 50 ==0){
			cout << "Event number: " << data_num << endl;
		}
		
        for ( int pix = 0; pix < 1440; pix++ ){
            hStartCell->Fill( pix, data_offset[pix] );
        }

        for ( unsigned int sl = 0; sl < data_roi; sl++){
            value = getValue(sl, pixelnr);
            //printf("value = %f\n", value);
            Ameas[ sl ] = value;
            h[tAmeas].SetBinContent(sl, value);
            
        }

        computeN1mean( data_roi );
        removeSpikes( data_roi );

        for ( unsigned int sl = 0; sl < data_roi; sl++){
            hPixelCellData.Fill(sl, Vcorr[ sl ]);
        }   

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

		vector<DiscOut> *my =  discriminator( Vslide, 3.5, 100 );
		 for (int p=0; p<my->size(); p++ ){
				sp->Fill(my->at(p).maxVal);
      }	
		delete my;


		if ( spikeDebug ){
			for ( unsigned int sl = 0; sl < data_roi; sl++){
				h[tVslide].SetBinContent( sl, Vslide[sl] );
				h[tVcfd].SetBinContent( sl, Vcfd[sl] );
				h[tVcfd2].SetBinContent( sl, Vcfd2[sl] );
			}   
        }
/*		
		vector<int> * zeros = zerosearch( Vcfd2, -1, 10, 20 );
		if (zeros->size() == 0 ){
			continue;
        }
		// check value of Vside at zero position	
		for ( int i=0; i<zeros->size(); i++){
			cout << zeros->at(i) << ":\t" << Vslide[ zeros->at(i) ]<<endl;
//			sp->Fill(Vslide[zeros->at(i)]);
		}
*/
        
		if ( spikeDebug ){

            CW->cd( tAmeas + 1);
            gPad->SetGrid();
            h[tAmeas].Draw();

            CW->cd( tN1mean + 1);
            gPad->SetGrid();
            h[tN1mean].Draw();

            CW->cd( tVcorr + 1);
            gPad->SetGrid();
            h[tVcorr].Draw();

           	// CW->cd( tVtest + 1);
		    // gPad->SetGrid();
   		    // h[tVtest].Draw();

       		cFilter->cd( Ntypes - tVslide );
            cFilter->cd(1);
	        gPad->SetGrid();
   		    h[tVslide].Draw();
        
       		cFilter->cd( Ntypes - tVcfd );
	        cFilter->cd(2);
            gPad->SetGrid();
   		    h[tVcfd].Draw();

			TLine zeroline(0, 0, 1024, 0);
   			zeroline.SetLineColor(kBlue);
			zeroline.Draw();

            cFilter->cd( Ntypes - tVcfd2 );
            cFilter->cd(3);
            gPad->SetGrid();
            h[tVcfd2].Draw();

            zeroline.Draw();
        
			CW->Update();
            cFilter->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") break;
        }


	}
    
    TCanvas * cSpektrum = new TCanvas();
    cSpektrum->cd();
    sp->Draw();
    TCanvas * cStartCell = new TCanvas();
    cStartCell->cd();
    hStartCell->Draw();
    hPixelCellData.Draw();

	delete cStartCell;
	
	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++ ) h[ tVcorr ].SetBinContent( i, Vcorr[i] );
}

void computeN1mean( int Samples ){

// compute the mean of the left and right neighbors of a channel

    for( int i = 0; i < Samples; 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.;
        }
        h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]);
    }
} // 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 * data_roi;
    slice_pt = pixel_pt + slice;
    drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi;
    cal_pt    = pixel_pt + drs_cal_offset;

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

    return( vcal );
}

void BookHistos( int Samples ){
// booking and parameter settings for all histos

    h = new TH1F[ Ntypes ];

    for ( int type = 0; type < Ntypes; type++){

        h[ type ].SetBins(Samples, 0, Samples);
        h[ type ].SetLineColor(1);
        h[ type ].SetLineWidth(2);

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

        // set Y axis paras
        h[ type ].GetYaxis()->SetLabelSize(0.1);
        h[ type ].GetYaxis()->SetTitleSize(0.1);
        h[ type ].GetYaxis()->SetTitleOffset(0.3);
        h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
    }
    CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
    CW->Divide(1, 3);
    cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
    cFilter->Divide(1, 3);

    hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);

}
int FOpenDataFile(fits &datafile){

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

    //Get the size of the data column
    data_roi = datafile.GetUInt("NROI");  // Value from header
    data_px = datafile.GetUInt("NPIX");
    data_n = datafile.GetN("Data");         //Size of column "Data" = #Pixel x ROI
    
    //Set the sizes of the data vectors
    data.resize(data_n,0);
    data_offset.resize(data_px,0);
    
    //Link the data to variables
    datafile.SetRefAddress("EventNum", data_num);
    datafile.SetVecAddress("Data", data);
    datafile.SetVecAddress("StartCellData", data_offset);
    datafile.GetRow(0);
    
    cout << "Opening data file successful..." << endl;

    return datafile.GetNumRows() ;
}