#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 "fits.h"

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

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

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

#define NPIX  1440
#define NCELL 1024

// data access and treatment
#define FAD_MAX_SAMPLES 1024

int NEvents;
vector<int16_t> Data;	// vector, which will be filled with raw data
vector<int16_t> StartCells;	// vector, which will be filled with DRS start positions
unsigned int EventID;	// index of the current event
UInt_t RegionOfInterest;	// Width of the Region, read out of the DRS
UInt_t NumberOfPixels;	// Total number of pixel, read out of the camera
size_t PXLxROI;                                        // Size of column "Data" = #Pixel x ROI

int NBSLeve = 1000;

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


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> SpikeEst(FAD_MAX_SAMPLES);  // corrected Values
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"

void computeSpikeEstimator( int );

// histograms
const int Ntypes = 7;
const unsigned int  // arranged by Dominik
	tAmeas 	= 0, 
	tSpikeEst = 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 = 1;


int fana2( 
       char *datafilename = "../raw/20110916_025.fits",
       const char *drsfilename = "../raw/20110916_024.drs.fits",
       int pixelnr = 0,
       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 = NULL;
       NEvents = openDataFits(
	datafilename,
	&datafile,
	Data,
	StartCells,
	EventID,
	RegionOfInterest,
	NumberOfPixels,
	PXLxROI );
-
       printf( "number of events in file: %d\n", NEvents );
       if (NEvents == 0){
           cout << "return code of openDataFits:" << datafilename<< endl;
           cout << "is zero -> aborting." << endl;
           return 1;
       }

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

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

// Book the histograms
    BookHistos( RegionOfInterest );

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

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

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

       datafile->GetRow( ev );
       if (ev % 50 ==0){
           cout << "Event number: " << EventID << endl;
       }

       // get the data of this pixel from the  Data  vector
       // apply the Drs Calibration and cut off 12 slices at the beginning
       // and at the end.
       calib_RC = applyDrsCalibration( Ameas, pixelnr,12,12, Offset, Gain, TriggerOffset,
	RegionOfInterest, Data, StartCells);
       if (calib_RC == 0){
           break;
       }
       
       computeSpikeEstimator( RegionOfInterest );
       removeSpikes( Ameas, Vcorr );
       sliding_avg( Vcorr, Vslide, 8 );

       //      for ( unsigned int sl = 0; sl < RegionOfInterest; 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);

       if ( spikeDebug ){
           for ( unsigned int sl = 0; sl < RegionOfInterest; sl++ ){

               h[tAmeas].SetBinContent( sl, Ameas[ sl ] );
               h[tVcorr].SetBinContent( sl, Vcorr[ sl ] );

               h[tVslide].SetBinContent( sl, Vslide[ sl ] );
               h[tVcfd].SetBinContent(   sl, Vcfd[ sl ] );
               h[tVcfd2].SetBinContent( sl, Vcfd2[ sl ] );
           }   
       }
		
       if ( spikeDebug ){
            CW->cd( tAmeas + 1);
            gPad->SetGrid();
            h[tAmeas].Draw();

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

            CW->cd( tVcorr + 1);
            gPad->SetGrid();
            h[tVcorr].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;
        }


     }
    
     return( 0 );
}

void computeSpikeEstimator( int Samples ){

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

    for( int i = 1; i < Samples-1; i++){

       N1mean[ i ] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
       SpikeEst[ i ]   = Ameas[ i ] - N1mean[ i ];
       h[tSpikeEst].SetBinContent( i, SpikeEst[ i ] );
    }
} // end of computeN1mean computation


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

}