/********************** 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 "openFits.h"
#include "openFits.c"

#include "discriminator.h"
#include "discriminator.C"
#include "zerosearch.h"
#include "zerosearch.C"
#include "factfir.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 TriggerOffsetROI, RC;
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

typedef struct{
  float maxAmpl;
  float countOfMax;
  } OverlayMaximum;
vector<OverlayMaximum> SingleMaximum;
// 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
TH1F *hMaxPeakOverlay;
//TH2F *hPixelPeakOverlay;
//TH2F *hEventPeakOverlay;
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/tpeak_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,310,400,300);
          cFiltered->Divide(1, 3);
        }
        // Canvases to show the peak template
        TCanvas *cPeakOverlay = NULL;
        cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 1, 1, 1200, 800);
        cPeakOverlay->Divide(2,1);
//        // All peaks of one Pixel
//        TCanvas *cPixelPeakOverlay = NULL;
//        cPixelPeakOverlay = new TCanvas("cPixelPeakOverlay", "Overlay of detected Peaks of all Events of one Pixel", 401, 1, 400, 400);
//        // All Peaks of one Event
//        TCanvas *cYProjection = NULL;
//        cYProjection = new TCanvas("cYProjection", "Y Projection of peak overlay for each bin ", 802, 1, 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 = openDataFits( 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
   RC = openCalibFits( drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, TriggerOffsetROI);
   if (RC == 0){
     cout << "return code of openCalibFits:" << drsfilename << endl;
     cout << "is zero -> aborting." << endl;
     return 1;
   }
// 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 , 8);
       // 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++){
           if (it->maxPos < 12 || it->maxPos > RegionOfInterest-12)
             continue;
           //domstest->Fill(it->maxPos);
           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), Vslide[pos]) ;
//              hPixelPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ;
//              hEventPeakOverlay->Fill( pos - (it->maxPos), 2*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(1);
                                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(1);
                    hPeakOverlay->Draw("COLZ");
                    cPeakOverlay->cd(2);
                  //  hMaxPeakOverlay->Draw("");
                    cPeakOverlay->Modified();
                    cPeakOverlay->Update();
                  }
				}
                if (breakout)
                        break;
        }

//-------------------------------------
// end of loop over Events
//-------------------------------------

//-------------------------------------
// Histogramm of Maximas in Overlay Spectra
//-------------------------------------
      // resize vector SingleMaximum to number of timeslices in Overlay Spectra
      SingleMaximum.resize(hPeakOverlay->GetNbinsX());

      //put maximumvalue of every Y-projection of every timeslice into vector
      for (int TimeSlice = 0; TimeSlice <= hPeakOverlay->GetNbinsX(); TimeSlice++ ){
        TH1D* hProjPeak =	hPeakOverlay->ProjectionY("proj_py", TimeSlice, TimeSlice, "");
        SingleMaximum[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 54;
        SingleMaximum[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() );
        hMaxPeakOverlay->SetBinContent(TimeSlice, SingleMaximum[ TimeSlice ].maxAmpl );
      }
      hMaxPeakOverlay->Fit("landau", "", "", -50, 250);

//-------------------------------------
// Draw Histograms
//-------------------------------------
        if (ProduceGraphic){
            //OVERLAY PEAKS
            cPeakOverlay->cd(1);
            hPeakOverlay->Draw("COLZ");
            cPeakOverlay->cd(2);
            hMaxPeakOverlay->Draw("");
            cPeakOverlay->Modified();
            cPeakOverlay->Update();

//            cPixelPeakOverlay->cd();
//            hPixelPeakOverlay->Draw("COLZ");
//            cPixelPeakOverlay->Modified();
//            cPixelPeakOverlay->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.)");
        }
/*
        ProjHistos = new TH1D[ NumberOfTimeSlices ];
        for ( int type = 0; type < NumberOfTimeSlices; type++){
                ProjHistos[ type ].SetBins(1024, 0, 1024);
                ProjHistos[ type ].SetLineColor(1);
                ProjHistos[ type ].SetLineWidth(2);

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

                // set Y axis paras
                ProjHistos[ type ].GetYaxis()->SetLabelSize(0.1);
                ProjHistos[ type ].GetYaxis()->SetTitleSize(0.1);
                ProjHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
                ProjHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
        }
*/
      hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks",
               hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
               512, -55.5, 200.5 );
        hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
        hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
        //hPeakOverlay->SetBit(TH2F::kCanRebin);
        hList.Add( hPeakOverlay );

      hMaxPeakOverlay = new TH1F("hMaxPeakOverlay", "Single Peak derived with Maximum of above Spektrum",
               hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
               (-1*hPeakOverlayXaxisLeft)-0.5,
               hPeakOverlayXaxisRight-0.5
                );
        hMaxPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
        hMaxPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
        hList.Add( hMaxPeakOverlay );

//        hPixelPeakOverlay = new TH2F("hPixelPeakOverlay", "Maximum of Statistic of overlayed Peaks",
//               hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
//               512, -55.5, 200.5 );
//        hPixelPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
//        hPixelPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
//        hList.Add( hPixelPeakOverlay );

//        hEventPeakOverlay = new TH2F("hEventPeakOverlay", "Overlay of detected Peaks of all Pixel of one Event",
//               hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
//               4096, -48.5, 200.5 );
//        hEventPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
//        hEventPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
//        hList.Add( hEventPeakOverlay );

}

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 )