/********************** FPulseTemplate ***********************
* read FACT raw data
* remove spikes
* calculate baseline
* find peaks (CFD and normal discriminator)
* compute pulse height and pulse integral spektrum of the peaks
+ search for Peaks in data
+ put peaks into different histos depending und Amplitude
+ draw histos for single, double, tripple, ....above 100mV Photonpulses
+ draw Parameterdevelopment depending on photon quantity
+ form a template (shape depending on number of photons)
  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 <TString.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 global 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;

//-----------------------------------------------------------------------------
// Initialisation of histograms
//-----------------------------------------------------------------------------

TH1F* h;
TH1F *debugHistos;
TH2F *hOverlayTemp = NULL;
TH2F *hSinglePeakOverlay;//histogrammm for overlay of detected Peaks
TH2F *hDoublePeakOverlay;
TH2F *hTripplePeakOverlay;
TH2F *hLargePeakOverlay;

TH1F *hMaximumTemp = NULL;
TH1F *hSinglePeakMaximum;
TH1F *hDoublePeakMaximum;
TH1F *hTripplePeakMaximum;
TH1F *hLargePeakMaximum;

int hPeakOverlayXaxisLeft,hPeakOverlayXaxisRight; 

TObjArray hList;
TObjArray hListBaseline;

//-----------------------------------------------------------------------------
// Functions
//-----------------------------------------------------------------------------

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

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
// MAIN - Funtion
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
int FPulseTemplate(
  char *datafilename    = "../data/2011/11/09/20111109_006.fits.gz",
  const char *drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz",
  const char *OutRootFileName = "../analysis/20111109_006.fits.gz.peaktypes.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 = 250,
  int verbosityLevel = 1, // different verbosity levels can be implemented here
  bool ProduceGraphic = true
  )

{
//-----------------------------------------------------------------------------
//	histogramm Settings
//-----------------------------------------------------------------------------

  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 *cPulsetypes = NULL;
		cPulsetypes = new TCanvas("cPulsetypes", "Overlay of detected Peaks", 1, 1, 1400, 700);
		cPulsetypes->Divide(4,2);

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

//-----------------------------------------------------------------------------
// Fill Overlay Histos
//-----------------------------------------------------------------------------
       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() ;
		   if (Vslide[it->maxPos] >= 5 && Vslide[it->maxPos] < 13) hOverlayTemp = &*hSinglePeakOverlay;
		   else if (Vslide[it->maxPos] >= 13 && Vslide[it->maxPos] < 23) hOverlayTemp = hDoublePeakOverlay;
		   else if (Vslide[it->maxPos] >= 23 && Vslide[it->maxPos] < 33) hOverlayTemp = hTripplePeakOverlay;
		   else if (Vslide[it->maxPos] >= 33) hOverlayTemp = hLargePeakOverlay;

           for ( int pos = Left; pos < Right; pos++){
//			 if (Vslide[it->maxPos] >= 5 && Vslide[it->maxPos] < 15) hSinglePeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]);
//			 else if (Vslide[it->maxPos] >= 15 && Vslide[it->maxPos] < 25) hDoublePeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]);
//			 else if (Vslide[it->maxPos] >= 25 && Vslide[it->maxPos] < 35) hTripplePeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]);
//			 else if (Vslide[it->maxPos] >= 35) hOverlayTemp = hLargePeakOverlay;
				hOverlayTemp->Fill( pos - (it->maxPos), Vslide[pos]) ;
           }
       }

//-----------------------------------------------------------------------------
// Spike Debug
//-----------------------------------------------------------------------------
       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
								cPulsetypes->cd(1);
								gPad->SetGrid();
								hSinglePeakOverlay->Draw("COLZ");

								cPulsetypes->cd(2);
								gPad->SetGrid();
								hDoublePeakOverlay->Draw("COLZ");

								cPulsetypes->cd(3);
								gPad->SetGrid();
								hTripplePeakOverlay->Draw("COLZ");

								cPulsetypes->cd(4);
								gPad->SetGrid();
								hLargePeakOverlay->Draw("COLZ");

								cPulsetypes->Modified();
								cPulsetypes->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 % 500 == 0){

                    //OVERLAY PEAKS
						  cPulsetypes->cd(1);
						  gPad->SetGrid();
						  hSinglePeakOverlay->Draw("COLZ");

						  cPulsetypes->cd(2);
						  gPad->SetGrid();
						  hDoublePeakOverlay->Draw("COLZ");

						  cPulsetypes->cd(3);
						  gPad->SetGrid();
						  hTripplePeakOverlay->Draw("COLZ");

						  cPulsetypes->cd(4);
						  gPad->SetGrid();
						  hLargePeakOverlay->Draw("COLZ");
				  //  cPulsetypes->cd(5);
				  //  hSinglePeakMaximum->Draw("");
					cPulsetypes->Modified();
					cPulsetypes->Update();
                  }
				}
                if (breakout)
                        break;
        }

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

//-------------------------------------
// Histogramm of Maximas in Overlay Spectra
//-------------------------------------
cout << "producing...Maximum peaks" << endl;
   for (unsigned int cnt = 1 ; cnt <= 4 ; cnt ++){
	 cout << "peak type: " << cnt << endl;
	 if (cnt == 1){
	   hMaximumTemp = hSinglePeakMaximum;
	   hOverlayTemp = hSinglePeakOverlay;
	 }
	 else if (cnt == 2){
	   hMaximumTemp = hDoublePeakMaximum;
	   hOverlayTemp = hDoublePeakOverlay;
	 }
	 else if (cnt == 3){
	   hMaximumTemp = hTripplePeakMaximum;
	   hOverlayTemp = hTripplePeakOverlay;
	 }
	 else if (cnt == 4){
	   hMaximumTemp = hLargePeakMaximum;
	   hOverlayTemp = hLargePeakOverlay;
	 }
   // resize vector SingleMaximum to number of timeslices in Overlay Spectra
	  SingleMaximum.resize(hOverlayTemp->GetNbinsX());
      //put maximumvalue of every Y-projection of every timeslice into vector
	  for (int TimeSlice = 0; TimeSlice <= hOverlayTemp->GetNbinsX(); TimeSlice++ ){
		TString histotitle;
		histotitle += "hproj_py";
		histotitle += cnt ;
		TH1D* hProjPeak =	hOverlayTemp->ProjectionY(histotitle, TimeSlice, TimeSlice, "");
		SingleMaximum[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 3.5;
        SingleMaximum[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() );
		hMaximumTemp->SetBinContent(TimeSlice, SingleMaximum[ TimeSlice ].maxAmpl );
      }
	  hMaximumTemp->Fit("landau", "", "", -50, 250);
	}
//-------------------------------------
// Draw Histograms
//-------------------------------------
        if (ProduceGraphic){
            //OVERLAY PEAKS
cout << "producing...So" << endl;
		  cPulsetypes->cd(1);
		  gPad->SetGrid();
		  hSinglePeakOverlay->ResetStats();
		  hSinglePeakOverlay->Draw("COLZ");
cout << "producing...DO" << endl;
		  cPulsetypes->cd(2);
		  gPad->SetGrid();
		  hDoublePeakOverlay->ResetStats();
		  hDoublePeakOverlay->Draw("COLZ");
cout << "producing...TO" << endl;
		  cPulsetypes->cd(3);
		  gPad->SetGrid();
		  hTripplePeakOverlay->ResetStats();
		  hTripplePeakOverlay->Draw("COLZ");
cout << "producing...LO" << endl;
		  cPulsetypes->cd(4);
		  gPad->SetGrid();
		  hLargePeakOverlay->ResetStats();
		  hLargePeakOverlay->Draw("COLZ");
cout << "producing...SM" << endl;
		  cPulsetypes->cd(5);
		  hSinglePeakMaximum->ResetStats();
		  hSinglePeakMaximum->Draw("");
cout << "producing...DM" << endl;
		  cPulsetypes->cd(6);
		  hDoublePeakMaximum->ResetStats();
		  hDoublePeakMaximum->Draw("");
cout << "producing...TM" << endl;
		  cPulsetypes->cd(7);
		  hTripplePeakMaximum->ResetStats();
		  hTripplePeakMaximum->Draw("");
cout << "producing...LM" << endl;
		  cPulsetypes->cd(8);
		  hLargePeakMaximum->ResetStats();
		  hLargePeakMaximum->Draw("");
cout << "producing...Done" << endl;




		  cPulsetypes->Modified();
		  cPulsetypes->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.)");
        }
*/
//-----------------------------------------------------------------------------
	  hSinglePeakOverlay = new TH2F("hSinglePeakOverlay", "Overlay of detected Single Photon Peaks",
               hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
			   512, -55.5, 300.5);
	  hSinglePeakOverlay->SetAxisRange(-5.5, 35.5, "Y");
		hSinglePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
		hSinglePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
		//hSinglePeakOverlay->SetBit(TH2F::kCanRebin);
		hList.Add( hSinglePeakOverlay );

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

		hDoublePeakOverlay = new TH2F("hDoublePeakOverlay", "Overlay of detected Double Photon Peaks",
		hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
		512, -55.5, 300.5);
		hDoublePeakOverlay->SetAxisRange(-5.5, 35.5, "Y");
		hDoublePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
		hDoublePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
		//hSinglePeakOverlay->SetBit(TH2F::kCanRebin);
		hList.Add( hDoublePeakOverlay );;

		hDoublePeakMaximum = new TH1F("hSinglePeakMaximum", "Double Peak derived with Maximum of above Spektrum",
				 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
				 (-1*hPeakOverlayXaxisLeft)-0.5,
				 hPeakOverlayXaxisRight-0.5
				  );
		hDoublePeakMaximum->SetAxisRange(-0.5, 25.5, "Y");
		  hDoublePeakMaximum->GetXaxis()->SetTitle( "Timeslices" );
		  hDoublePeakMaximum->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
		  hList.Add( hDoublePeakMaximum );
  //-----------------------------------------------------------------------------

		hTripplePeakOverlay= new TH2F("hTripplePeakOverlay", "Overlay of detected Tripple Photon Peaks",
		hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
		512, -55.5, 300.5);
		hTripplePeakOverlay->SetAxisRange(-5.5, 35.5, "Y");
		hTripplePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
		hTripplePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
		//hTripplePeakOverlay->SetBit(TH2F::kCanRebin);
		hList.Add( hTripplePeakOverlay );;

		hTripplePeakMaximum = new TH1F("hSinglePeakMaximum", "Triple Peak derived with Maximum of above Spektrum",
				 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
				 (-1*hPeakOverlayXaxisLeft)-0.5,
				 hPeakOverlayXaxisRight-0.5
				  );
		hTripplePeakMaximum->SetAxisRange(-0.5, 25.5, "Y");
		  hTripplePeakMaximum->GetXaxis()->SetTitle( "Timeslices" );
		  hTripplePeakMaximum->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
		  hList.Add( hTripplePeakMaximum );
  //-----------------------------------------------------------------------------

		hLargePeakOverlay= new TH2F("hLargePeakOverlay", "Overlay of detected Multi Photon Peaks",
		hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
		512, -5.5, 300.5 );
		hLargePeakOverlay->SetAxisRange(-5.5, 200.5, "Y");
		hLargePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
		hLargePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
		//hLargePeakOverlay->SetBit(TH2F::kCanRebin);
		hList.Add( hLargePeakOverlay );;

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

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