/********************** 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;
bool UseThisPeak=false;
int NEvents;
vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;
unsigned int CurrentEventID;
size_t PXLxROI;
UInt_t RegionOfInterest;
UInt_t NumberOfPixels;
TString histotitle;

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;

typedef struct{
  TString name;
  TString title;
  TString Mname;
  TString Mtitle;
  float yMax;
} hist_t;

hist_t histSetting[4];

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= NULL;
TH1F *debugHistos= NULL;
TH2F *hOverlayTemp = NULL;
TH2F *hPeakOverlay[4];//histogrammm for overlay of detected Peaks
TH2F *hSinglePeakOverlay;//histogrammm for overlay of detected Peaks
TH2F *hDoublePeakOverlay;
TH2F *hTripplePeakOverlay;
TH2F *hLargePeakOverlay;

TH1F *hMaximumTemp = NULL;
TH1F *hPeakMaximum[4];
TH1F *hSinglePeakMaximum;
TH1F *hDoublePeakMaximum;
TH1F *hTripplePeakMaximum;
TH1F *hLargePeakMaximum;
TH1D *hProjPeak = NULL;
int gPeakOverlayXaxisLeft;
int hPeakOverlayXaxisRight;

TObjArray hList;
TObjArray hListBaseline;

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

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

//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
// MAIN - Funtion
//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
int FPulseTemplate(
  char*         datafilename    = "../data/2011/11/10/20111110_005.fits.gz",
  const char*   drsfilename     = "../data/2011/11/10/20111110_003.drs.fits.gz",
  const char*   OutRootFileName = "../analysis/20111110_003.test.root",
  bool          ProduceGraphic  = true,
  bool          spikeDebug      = false,
  bool          fitdata         = false,
  int           verbosityLevel  = 1, // different verbosity levels can be implemented here
  int           firstevent      = 0,
  int           nevents         = -1,
  int           firstpixel      = 400,
  int           npixel          = 1,
  int           PintervalWidth  = 5,
  int           PintervalMin    = 8,
  int           PintervalMax    = 32, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
  int           avg1            = 14,
  int           avg2            = 8,
  int           OverlayWindowLeft = 70,
  int           OverlayWindowRight = 230)
{
//-----------------------------------------------------------------------------
//	histogramm Settings
//-----------------------------------------------------------------------------

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


        for (int stepper = 0; stepper < 4; stepper++){
          if (stepper == 0) histSetting[stepper].yMax = PintervalMin;
          else if (stepper == 3) histSetting[stepper].yMax = PintervalMax;
          else histSetting[stepper].yMax = ((PintervalMax-PintervalMin)*(stepper)/3)+PintervalMin;
          cout << "Max @ " << histSetting[stepper].yMax << endl;
          histSetting[stepper].name = "hPeakOverlay";
          histSetting[stepper].name += stepper;
          histSetting[stepper].title = "PeakOverlay with Amplitude around ";
          histSetting[stepper].title += histSetting[stepper].yMax;
          histSetting[stepper].title += " mV";
          histSetting[stepper].Mname = "hPeakMaximum";
          histSetting[stepper].Mname += stepper;
          histSetting[stepper].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of ";
          histSetting[stepper].Mtitle += histSetting[stepper].yMax;
          histSetting[stepper].Mtitle += " mV";
        }

//-----------------------------------------------------------------------------
// 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){
             cout << endl;
             if (pix == firstpixel){
                 cout << "...processing Event: " << CurrentEventID << "/" << nevents << endl;
             }
           }

//-------------------------------------
// Apply Filters
//-------------------------------------
       cout << "...applying DrsCalibration";
       applyDrsCalibration( Ameas,pix,12,12,
           drs_basemean, drs_gainmean, drs_triggeroffsetmean,
           RegionOfInterest, AllPixelDataVector, StartCellVector);
       cout << "...done " << endl;

       // finds spikes in the raw data, and interpolates the value
       // spikes are: 1 or 2 slice wide, positive non physical artifacts
       cout << "...removeing Spikes";
       removeSpikes (Ameas, Vcorr);
       cout << "...done " << endl;
       // filter Vcorr with sliding average using FIR filter function
       cout << "...applying sliding average filter";
       sliding_avg(Vcorr, Vslide, avg1);
       cout << "...done " << endl;
       // filter Vslide with CFD using FIR filter function
       cout << "...apllying factfir filter";
       factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
       cout << "...done " << endl;
       // filter Vcfd with sliding average using FIR filter function
       cout << "...applying 2nd sliding average filter";
       sliding_avg(Vcfd, Vcfd2, avg2);
       cout << "...done " << endl;

//-------------------------------------
// Search vor Zero crossings
//-------------------------------------
       // 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() ;

           cout << "...choosing Histogram" ;
           for(int stepper = 0; stepper < 4; stepper++){
             if (Vslide[it->maxPos] >= (histSetting[stepper].yMax - (PintervalWidth/2)) && Vslide[it->maxPos] < (histSetting[stepper].yMax + PintervalWidth/2)){
               hOverlayTemp = hPeakOverlay[stepper];
               cout << "...#" << stepper ;
               UseThisPeak = true;
               break;
             }
             else if (stepper > 2){
               UseThisPeak = false;
               cout << "...NONE" << endl ;
             }
           }

//		   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;
            if (UseThisPeak){
            cout << "...filling" ;
             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]) ;
             }
             cout << "...done" << endl;
            }

       }

//-----------------------------------------------------------------------------
// 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
                                for(int stepper = 0; stepper < 4; stepper++){
                                cPulsetypes->cd(stepper+1);
                                gPad->SetGrid();
                                hPeakOverlay[stepper]->Draw("COLZ");
                                }


//								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 % 50 == 0){

                    //OVERLAY PEAKS
                          cout << "...drawing histograms after Loop over Pixels" ;
                          for(int stepper = 0; stepper < 4; stepper++){
                          cPulsetypes->cd(stepper+1);
                          gPad->SetGrid();
                          hPeakOverlay[stepper]->Draw("COLZ");
                          cout << "..." << stepper << "/3" ;
                          }
//						  cPulsetypes->cd(1);
//						  gPad->SetGrid();
//						  hSinglePeakOverlay->Draw("COLZ");
//
//						  cPulsetypes->cd(2);
//						  gPad->SetGrid();
//						  hDoublePeakOverlay->Draw("hLargePeakOverlayCOLZ");
//
//						  cPulsetypes->cd(3);
//						  gPad->SetGrid();
//						  hTripplePeakOverlay->Draw("COLZ");
//
//						  cPulsetypes->cd(4);
//						  gPad->SetGrid();
//						  hLargePeakOverlay->Draw("COLZ");
                  //  cPulsetypes->cd(5);
                  //  hSinglePeakMaximum->Draw("");
                          cout << "...done" << endl;
                    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 = 0 ; cnt < 4 ; cnt ++){
     cout << "\t ...peak type " << cnt;
     hOverlayTemp = hPeakOverlay[cnt];
     hMaximumTemp = hPeakMaximum[cnt];
//	 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());
      cout << "\t ...peak type " << cnt;
      //put maximumvalue of every Y-projection of every timeslice into vector
      for (int TimeSlice = 0; TimeSlice <= hOverlayTemp->GetNbinsX(); TimeSlice++ ){
        histotitle = "hproj_py";
        histotitle += cnt ;
        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 );
      }
      cout << "...done" << endl;
      if (fitdata){
        cout << "...calculating Landaufit" ;
        hMaximumTemp->Fit("landau", "", "", -50, 250);
        cout << "...done" << endl;
        }
    }
//-------------------------------------
// Draw Histograms
//-------------------------------------
        if (ProduceGraphic){
            //OVERLAY PEAKS
          cout << "...drawing overlay histograms" ;
          for(int stepper = 0; stepper < 4; stepper++){
          cPulsetypes->cd(stepper+1);
          gPad->SetGrid();
          hPeakOverlay[stepper]->ResetStats();
          hPeakOverlay[stepper]->Draw("COLZ");
          }
          cout << "...done" << endl;
          cout << "...drawing peak-shape histograms" ;
          for(int stepper = 0; stepper < 4; stepper++){
          cPulsetypes->cd(5+stepper);
          gPad->SetGrid();
          hPeakMaximum[stepper]->ResetStats();
          hPeakMaximum[stepper]->Draw();
          }
          cout << "...done" << endl;

//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.)");
        }
*/

//-----------------------------------------------------------------------------

        for (int stepper = 0; stepper < 4; stepper ++){
      hPeakOverlay[stepper] = new TH2F(histSetting[stepper].name, histSetting[stepper].title,
               gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
               512, -55.5, 300.5);
      hPeakOverlay[stepper]->SetAxisRange(-5.5, histSetting[stepper].yMax+ 25.5, "Y");
        hPeakOverlay[stepper]->GetXaxis()->SetTitle( "Timeslices" );
        hPeakOverlay[stepper]->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
        //hSinglePeakOverlay->SetBit(TH2F::kCanRebin);
        hList.Add( hPeakOverlay[stepper] );

      hPeakMaximum[stepper] = new TH1F(histSetting[stepper].Mname, histSetting[stepper].Mtitle,
               gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
               (-1*gPeakOverlayXaxisLeft)-0.5,
               hPeakOverlayXaxisRight-0.5
                );
      hPeakMaximum[stepper]->SetAxisRange(-0.5, histSetting[stepper].yMax + 15.5, "Y");
        hPeakMaximum[stepper]->GetXaxis()->SetTitle( "Timeslices" );
        hPeakMaximum[stepper]->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
        hList.Add( hPeakMaximum[stepper] );
      }
//-----------------------------------------------------------------------------



//-----------------------------------------------------------------------------
      hSinglePeakOverlay = new TH2F("hSinglePeakOverlay", "Overlay of detected Single Photon Peaks",
               gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-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",
               gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
               (-1*gPeakOverlayXaxisLeft)-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",
        gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-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",
                 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
                 (-1*gPeakOverlayXaxisLeft)-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",
        gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-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",
                 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
                 (-1*gPeakOverlayXaxisLeft)-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",
        gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-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",
                 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
                 (-1*gPeakOverlayXaxisLeft)-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",
//               gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-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",
//               gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-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 )