/********************** Peak-Template ***********************

this makro shell search for Peaks in data and form a template
of these, so it can be used for detection of other peaks

*************************************************************/

//root libraries

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

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

#define HAVE_ZLIB

//rootmacros

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

//-----------------------------------------------------------------------------
// Decleration of variables
//-----------------------------------------------------------------------------
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 drs_n;
vector<float> drs_basemean;
vector<float> drs_gainmean;
vector<float> drs_triggeroffsetmean;

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 *hPeakOverlay;  //histogrammm for overlay of detected Peaks
int hPeakOverlayXaxisLeft,hPeakOverlayXaxisRight; 

TObjArray hList;
TObjArray hListBaseline;

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

//-----------------------------------------------------------------------------
// Main
//-----------------------------------------------------------------------------
int tpeak(
  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    = 14,
  int avg2    = 8,
	int OverlayWindowLeft = 50,
	int OverlayWindowRight = 150,
  int verbosityLevel = 1, // different verbosity levels can be implemented here
  bool ProduceGraphic = true
  )

{
hPeakOverlayXaxisLeft = OverlayWindowLeft;
hPeakOverlayXaxisRight = OverlayWindowRight;

	gStyle->SetPalette(1,0);
	gROOT->SetStyle("Plain");

//-----------------------------------------------------------------------------
// Create (pointer to) Canvases, which are used in every run,
// also in 'non-debug' runs
//-----------------------------------------------------------------------------
        // Canvases only need if spike Debug, but I want to deklare
        // the pointers anyway ...
        TCanvas *cFiltered = NULL;
        if (spikeDebug){
                cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",1,10,400,300);
                cFiltered->Divide(1, 3);
        }
        // Canvases to show the peak template
        TCanvas *cPeakOverlay = NULL;
        cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 1, 310, 400, 400);


//-----------------------------------------------------------------------------
// read FACT raw data
// 	* remove spikes
//	* calculate baseline
//	* find peaks (CFD and normal discriminator)
//	* compute pulse height and pulse integral spektrum of the peaks
//-----------------------------------------------------------------------------


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

//-----------------------------------------------------------------------------
// prepare datafile
//-----------------------------------------------------------------------------

// Open the data file

  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 = OpenDataFile( 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;

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

// Book the histograms
        BookHistos( );

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

        //-------------------------------------
        // Loop over every Pixel of Event
        //-------------------------------------

                for ( int pix = firstpixel; pix < firstpixel + npixel; pix++ ){

												if (verbosityLevel > 0){
									        if (pix == firstpixel){
									          cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
									        }
									      }

											applyDrsCalibration( Ameas,pix,12,12,
								        drs_basemean, drs_gainmean, drs_triggeroffsetmean,
								        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 ... 10 is okay as well
                        vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
                        // zXings means "zero cross ings"
                        EnlargeRegion(*zXings, 10, 10);
                        findAbsMaxInRegions(*zXings, Vslide);
									      removeMaximaBelow( *zXings, 3.0);
									      removeRegionWithMaxOnEdge( *zXings, 2);
									      removeRegionOnFallingEdge( *zXings, 100);

                        //following Code should produce the Overlay of peaks
                        vector<Region>::iterator it;

                                for (it = zXings->begin() ; it < zXings->end() ; it++){
                                  int Left = it->maxPos - OverlayWindowLeft;
                                  int Right = it->maxPos + OverlayWindowRight;
                                  if (Left < 0)
                                    Left =0;
                                  if (Right > (int)Vcorr.size() )
                                    Right=Vcorr.size() ;
                                        for ( int pos = Left; pos < Right; pos++){

                                            hPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ;
                                        }
                                }


                        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[Vcorr_].Draw();

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

                                cFiltered->Update();

                                //OVERLAY PEAKS
                                cPeakOverlay->cd();
                                hPeakOverlay->Draw("COLZ");
                                cPeakOverlay->Modified();
                                cPeakOverlay->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;
                        if (breakout)
                                break;

                }
        //-------------------------------------
        // end of loop over pixels
        //-------------------------------------
				if (ProduceGraphic){
						if (ev % 50 == 0){
                //OVERLAY PEAKS
                cPeakOverlay->cd();
                hPeakOverlay->Draw("COLZ");
                cPeakOverlay->Modified();
                cPeakOverlay->Update();
						}
				}
                if (breakout)
                        break;
        }
//-------------------------------------
// end of loop over Events
//-------------------------------------
				if (ProduceGraphic){
                //OVERLAY PEAKS
                cPeakOverlay->cd();
                hPeakOverlay->Draw("COLZ");
                cPeakOverlay->Modified();
                cPeakOverlay->Update();
				}

        SaveHistograms( OutRootFileName );

        return( 0 );
}

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

        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.)");
        }
        hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks", 
					hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
					4096, -48.5, 1999.5 );
        hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
        hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
        hList.Add( hPeakOverlay );

}

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
} // end of SaveHistograms( char * loc_fname )