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

#include <stdio.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"

//#include "zerosearch.C"
#include "factfir.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> 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



// 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 hPixelCellData(
    "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
TH1F *hMeanBsl, *hpltMeanBsl;
TH1F *hRmsBsl, *hpltRmsBsl; 
TObjArray hList;
TObjArray hListBaseline;

void BookHistos( int , int);
void SaveHistograms( const char * );

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


int fbsl( 
	const char *datafilename 		= "path-to-datafile.fits.gz",
	const char *drsfilename = "path-to-calibfile.drs.fits.gz",
	int firstevent 			= 0, 
	int nevents 			= -1, 
	int firstpixel		= 0,
	int npixel				= -1,
	const char *TextOutFileName = "",
	const char *RootOutFileName = "",
	bool produceGraphic = true
){
	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;
	if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels;

	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( RegionOfInterest, npixel );
	size_t calib_RC = 1;
 
		// loop over events
    for ( int ev = firstevent; ev < firstevent + nevents; ev++) {

	    datafile->GetRow( ev );
		
		if( ev % 100 == 0){
			cout << "Event ID: " << EventID << endl;
			}
	
		// loop over pixel
		for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){


			// 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,pix,12,12,
				Offset, Gain, TriggerOffset,
				RegionOfInterest, Data, StartCells);
			if (calib_RC == 0){
				break;
			}

			// 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
			// 8 is here the HalfWidth of the sliding average window.
			sliding_avg(Vcorr, Vslide, 8);

            for ( unsigned int sl = 0; sl < Vslide.size() ; sl++){
               // hPixelCellData.Fill( sl, Vcorr[sl] );
               hBaseline[pix-firstpixel]->Fill( Vslide[sl] ) ;
            }   
	    }

			if (calib_RC == 0){
				break;
			}
	
	}


			if (calib_RC == 0){
				cout << "DRS Calibration didn't work ... aborting" << endl;
				return 2;
			}

    FILE *fp;
    TString fName; 

    fName = TextOutFileName;
	if (fName.Length() != 0){
		 cout << "saving text to " << TextOutFileName << endl;

	
    fp = fopen(fName, "a+");
    fprintf( fp, "FILE: %s\n", datafilename );
    fprintf( fp, "NEVENTS: %d\n", nevents);
    NBSLeve = nevents; // this has to be changed when computation on a subset of a run is implemented
    fprintf( fp, "NBSLEVE:  %d\n", NBSLeve );

    for (int pix = firstpixel; pix < firstpixel+npixel; pix++) {
        float vmaxprob = hBaseline[pix-firstpixel]->GetXaxis()->GetBinCenter(
            hBaseline[pix-firstpixel]->GetMaximumBin() );

        fprintf( fp, "%8.3f", vmaxprob );

        hMeanBsl->Fill( vmaxprob );
        hpltMeanBsl->SetBinContent(pix+1, vmaxprob );

        hRmsBsl->Fill(hBaseline[pix-firstpixel]->GetRMS() );
        hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );   
    }
    fprintf( fp, "\n" );

    fclose( fp );
	}
	if (RootOutFileName != ""){
		cout << "saving histograms to " << RootOutFileName << endl;
    SaveHistograms( RootOutFileName );
	}
		if (produceGraphic){
    	TCanvas * cMeanBsl = new TCanvas();
    	cMeanBsl->cd();
    	hMeanBsl->Draw();
    	cMeanBsl->Update();

    	TCanvas * cRmsBsl = new TCanvas();
    	cRmsBsl->cd();
    	hRmsBsl->Draw();
    	cMeanBsl->Update();
		}
	return( 0 );
}

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

    // histograms for baseline extraction
    char hName[500];
    char hTitle[500];

    TH1F *h;
    
    printf("inside BookHistos\n");

    for( int i = 0; i < NumberOfPixel; i++ ) {

        // printf("call sprintf [%d]\n", i );
        sprintf(&hTitle[0],"all events all slices of pixel %d", i);
        sprintf(&hName[0],"base%d", i);
        // printf("call sprintf [%d] done\n", i );

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

        // printf("histo booked\n");
        h->GetXaxis()->SetTitle( "Sample value (mV)" );
        h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
        // printf("histo title set\n");
        hListBaseline.Add( h );
        // printf("histo added to List\n");
        hBaseline[i] = h;
        // printf("histo assigned to array\n");
    }

    printf("made HBaseline * 1440\n");

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


void SaveHistograms( const 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
  //fName.Remove(fName.Length() - 5); // remove the extension .root
  //fName += "_histo.root"; // add the new extension
  //fName += ".root"; 

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

  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.Close(); // close the file

} // end of function: void ana::SaveHistograms( void )