/********************** 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 <TProfile.h>

#include <stdio.h>
#include <stdint.h>
#include <cstdio>
#include <iostream>
#include <fstream>
using namespace std;

#define NPIX  1440
#define NCELL 1024
#define FAD_MAX_SAMPLES 1024
#define MAX_PULS_ORDER 1
#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 gNEvents;
vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;
unsigned int CurrentEventID;
size_t PXLxROI;
UInt_t gRegionOfInterest;
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

float gGainMean = 9;
float gBSLMean = 0;

typedef struct{
  double maxAmpl;
  int countOfMax;
  } OverlayMaximum;

//typedef struct{
//  TString name;
//  TString title;
//  TString xTitle;
//  TString yTitle;
//  float yMax;
//  float yMin;
//} histAttrib_t;

//histAttrib_t* gHistoAttrib[MAX_PULS_ORDER];
//histAttrib_t* gMaximumAttrib[MAX_PULS_ORDER];
//histAttrib_t* gMaxGausAttrib[MAX_PULS_ORDER];
//histAttrib_t* gProfileAttrib[MAX_PULS_ORDER];


// histograms
const int Number_Of_Debug_Histo_Types = 7;

const unsigned int
        Ameas_ 	= 0,
        N1mean_ = 1,
        Vcorr_  = 2,
        Vtest_  = 3,
        Vslide_ = 4,
        Vcfd_   = 5,
        Vcfd2_  = 6;

//const char* gHistoTypes[8] = {
//    "hPixelOverlay",
//    "hPixelMax",
//    "hPixelEdgeOverlay",
//    "hPixelProfile",
//    "hAllPixelOverlay",
//    "hAllPixelMax",
//    "hAllPixelMaxGaus",
//    "hAllPixelProfile"
//};

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

// Temporary Objects
TH1F*       debugHistos = NULL;
TH2F*       hOverlayTemp = NULL;
TH1F*       hMaximumTemp = NULL;
TH1D*       hProjPeak = NULL;
TH1F*       hTesthisto = NULL;
TH2F*       hTesthisto2 = NULL;

//Pixelwise Histograms
TH2F*       hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks
TH2F*       hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
TProfile*   hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
TProfile*   hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
TH1F*       hPixelMax[MAX_PULS_ORDER] = {NULL};
TH1F*       hPixelMedian[MAX_PULS_ORDER] = {NULL};
TH1F*       hPixelMean[MAX_PULS_ORDER] = {NULL};

//All Pixel Histograms
TH2F*       hAllPixelOverlay[MAX_PULS_ORDER] = {NULL};
TProfile*   hAllPixelProfile[MAX_PULS_ORDER] = {NULL};
TProfile*   hAllPixelProfile2[MAX_PULS_ORDER] = {NULL};
TH1F*       hAllPixelMax[MAX_PULS_ORDER]     = {NULL};
TH1F*       hAllPixelMedian[MAX_PULS_ORDER] = {NULL};
TH1F*       hAllPixelMean[MAX_PULS_ORDER] = {NULL};
TH2F*       hAllPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};

//Histogram Parameters
Int_t       gPixelOverlayXaxisLeft = 0;
Int_t       gPixelOverlayXaxisRight = 0;

//Root-File Objects
TObjArray*  hList = NULL;
TObjArray*  hAllPixelList = NULL;
TObjArray*  hRootList = NULL;

//Canvases
TCanvas*    cgpPixelPulses[MAX_PULS_ORDER] = {NULL};
TCanvas*    cgpAllPixelPulses[MAX_PULS_ORDER] = {NULL};
TCanvas*    cgpTestHistos = NULL;
//TCanvas*    gpcDevelopment = NULL;

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

void BookHistos(int );
void BookPixelHistos(int );
void BookTestHistos( int );

void DeletePixelHistos(int );
void SaveHistograms( const char*, const char*, TObjArray*, int );

void FillHistograms(vector<Region>*, int, int, int);
void DrawPulseHistograms(int, int);
void DrawMaximumHistograms( int, int );
void DrawTestHistograms( int);

void PlotPulseShapeOfMaxPropability(unsigned int, int, int);
void FitMaxPropabilityPuls( TH1F* , int );

Double_t MedianOfH1( TH1 *h1 );
void PlotMedianEachSliceOfPulse(unsigned int, int, int);

void UpdateCanvases( int, int);
void DeletePixelCanvases( int );

void CreateRootFile(const char*, int );
TFile *OpenRootFile(const char*, int );
void CloseRootFile(TFile *tf);

TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file
TString CreateSubDirName(const char* );

void WritePixelTemplateToCsv( TString, const char*, const char*, int, int);
void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int);
//void SetHistogrammSettings( const char*, int, int , int);
//void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int);
//void DeletePulseCanvas( TCanvas* , unsigned int);


//----------------------------------------------------------------------------
//----------------------------------------------------------------------------
// 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/FPulseTemplate/20111109_006/20111109_006.pulses.root",
  const char*   OutPutPath          = "../analysis/FPulseTemplate/20111109_006/Templates/",
  bool          ProduceGraphic      = false,
  int           refresh_rate        = 500,      //refresh rate for canvases
  bool          spikeDebug          = false,
  bool          debugPixel          = false,
  bool          fitdata             = false,
  int           verbosityLevel      = 0,        // different verbosity levels can be implemented here
  int           firstevent          = 0,
  int           nevents             = -1,
  int           firstpixel          = 0,
  int           npixel              = -1,
  int           AmplWindowWidth     = 14,       //Width of Window for selection of pulses to histograms
  float         GainMean           = 8,
  float         BSLMean            = -1,        //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
  int           avg1                = 8,
  int           avg2                = 8,
  int           OverlayWindowLeft   = 70,
  int           OverlayWindowRight  = 230
        )
{
    gGainMean = GainMean;
    gBSLMean = BSLMean;
//----------------------------------------------------------------------------
//	Save-Root-File Settings
//----------------------------------------------------------------------------
    CreateRootFile( OutRootFileName, verbosityLevel );


//----------------------------------------------------------------------------
//	root global Settings
//----------------------------------------------------------------------------

    gPixelOverlayXaxisLeft = OverlayWindowLeft;
    gPixelOverlayXaxisRight = OverlayWindowRight;

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

        for (
             int pulse_order = MAX_PULS_ORDER - 1;
             pulse_order >= 0 ;
             pulse_order--
              )
        {
            TString cName   ="cgpPixelPulses";
                    cName   += pulse_order;
            TString cTitle   ="Pulses of Order: ";
                    cTitle   += pulse_order;
            cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800);
            cgpPixelPulses[pulse_order]->Divide(2, 2);
                    cName   = "cgpAllPixelPulses";
                    cName   += pulse_order;
                    cTitle   ="AllPixel average of puls shapes of Order: ";
                    cTitle   += pulse_order;
            cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 );
            cgpAllPixelPulses[pulse_order]->Divide(2, 2);
        }

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

        cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 );
        cgpTestHistos->Divide(2,0);
//-----------------------------------------------------------------------------
// 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.
    gNEvents = openDataFits(
                    datafilename,
                    &datafile,
                    AllPixelDataVector,
                    StartCellVector,
                    CurrentEventID,
                    gRegionOfInterest,
                    NumberOfPixels,
                    PXLxROI,
                    verbosityLevel
                );

    if (gNEvents == 0)
    {
        cout << "return code of OpenDataFile:" << datafilename<< endl;
        cout << "is zero -> aborting." << endl;
        return 1;
    }

    if (verbosityLevel > 0)
    {
        cout << endl <<"number of events in file: "<< gNEvents << "\t";
    }

    if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all!

    if (verbosityLevel > 0)
    {
        cout <<"of, which "<< nevents << " will be processed"<< endl;
        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(verbosityLevel );
    BookTestHistos( verbosityLevel );
//-------------------------------------
// Loop over every Pixel
//-------------------------------------
    for ( int pixel = firstpixel; pixel < firstpixel + npixel; pixel++ )
    {
        if (verbosityLevel >= 0)
        {
            cout << "------------------------------------------------" << endl
                 << "...processing Pixel: " << pixel << endl;
        }

        BookPixelHistos(verbosityLevel );

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

            if (verbosityLevel > 0)
           {
             cout << "-------------------------------------" << endl
                  << "...processing Event: " << CurrentEventID
                  << "/" << nevents << " of Pixel: " << pixel << endl;
           }

//-------------------------------------
// Apply Calibration
//-------------------------------------
            if (verbosityLevel > 2) cout << "...applying DrsCalibration";
            applyDrsCalibration(
                        Ameas,
                        pixel,
                        12,
                        12,
                        drs_basemean,
                        drs_gainmean,
                        drs_triggeroffsetmean,
                        gRegionOfInterest,
                        AllPixelDataVector,
                        StartCellVector
                        );
            if (verbosityLevel > 2) cout << "...done " << endl;

//-------------------------------------
// Apply Filters
//-------------------------------------
            // finds spikes in the raw data, and interpolates the value
            // spikes are: 1 or 2 slice wide, positive non physical artifacts
            if (verbosityLevel > 2) cout << "...removeing Spikes";
            removeSpikes (Ameas, Vcorr);
            if (verbosityLevel > 2) cout << "...done " << endl;

            // filter Vcorr with sliding average using FIR filter function
            if (verbosityLevel > 2) cout << "...applying sliding average filter";
            sliding_avg(Vcorr, Vslide, avg1);
            if (verbosityLevel > 2) cout << "...done " << endl;

            // filter Vslide with CFD using FIR filter function
            if (verbosityLevel > 2) cout << "...apllying factfir filter";
            factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
            if (verbosityLevel > 2) cout << "...done " << endl;

            // filter Vcfd with sliding average using FIR filter function
            if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter";
            sliding_avg(Vcfd, Vcfd2, avg2);
            if (verbosityLevel > 2) cout << "...done " << endl;

//-------------------------------------
// Search vor Zero crossings
//-------------------------------------
            if (verbosityLevel > 2) cout << endl << "...searching 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>* pZXings = zerosearch( Vcfd2 , 1 , 8);
            // pZXings means "zero cross ings"
            EnlargeRegion(*pZXings, 10, 10);
            findAbsMaxInRegions(*pZXings, Vslide);
            removeMaximaBelow( *pZXings, 3.0);
            removeRegionWithMaxOnEdge( *pZXings, 2);
            removeRegionOnFallingEdge( *pZXings, 100);
            findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel );
            if (verbosityLevel > 2) cout << "...done" << endl;

//-----------------------------------------------------------------------------
// Fill Overlay Histos
//-----------------------------------------------------------------------------
        FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel );

//-----------------------------------------------------------------------------
// 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 < gRegionOfInterest; 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[Ameas_].Draw();

            cFiltered->cd(2);
            gPad->SetGrid();
            debugHistos[Vcorr_].Draw();

            cFiltered->cd(3);
            gPad->SetGrid();
            debugHistos[Vslide_].Draw();

            TBox *OneBox;
            vector<TBox*> MyBoxes;
            for (unsigned int i=0; i<pZXings->size(); i++){
                    OneBox = new TBox(
                            pZXings->at(i).maxPos -10 ,
                            pZXings->at(i).maxVal -0.5,
                            pZXings->at(i).maxPos +10 ,
                            pZXings->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();


            //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 pZXings;
        if (breakout)   break;
//-------------------------------------
// Draw 1. Set of Pixel Histograms
//-------------------------------------
            if ((ev % refresh_rate) == 0)
            {
                DrawPulseHistograms(
                        verbosityLevel,
                        MAX_PULS_ORDER
                        );
                DrawTestHistograms(
                        verbosityLevel
                        );
              if (ProduceGraphic)
              {
                UpdateCanvases(
                                verbosityLevel,
                                MAX_PULS_ORDER
                                );
              }
            }

        if (breakout) break;

        if (verbosityLevel > 2)
        {
            cout << endl << "-------------------------------------"<< endl;
        }
        }//End of Loop over Events
//-------------------------------------
// end of Loops over Events
//-------------------------------------




//-------------------------------------
// Histogramms of Maximas in Overlay Spectra
//-------------------------------------

         PlotPulseShapeOfMaxPropability(
                    MAX_PULS_ORDER,
                    fitdata,
                    verbosityLevel
                 );

         PlotMedianEachSliceOfPulse(
                     MAX_PULS_ORDER,
                     fitdata,
                     verbosityLevel
                 );

         WritePixelTemplateToCsv(
                     OutPutPath,
                     "PulseTemplate_PointSet",
                     "Maximum",
                     pixel,
                     verbosityLevel
                     );

       if (verbosityLevel > 2) cout << "...done" << endl;
       //here is what happends at the end of each loop over all pixels
       //Save Histograms of Pixel into Output rootfile
       DrawPulseHistograms(
               verbosityLevel,
               MAX_PULS_ORDER
               );

       DrawMaximumHistograms(
               verbosityLevel,
               MAX_PULS_ORDER
               );

        if (ProduceGraphic)
        {
        UpdateCanvases(
                       verbosityLevel,
                       MAX_PULS_ORDER
                       );
        }

//       SaveHistograms(
//                   OutRootFileName,
//                   CreateSubDirName(pixel),
//                   hList,
//                   verbosityLevel
//                   );

       if (debugPixel)
       {
           //Process gui events asynchronously during input
           cout << endl;
           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") {
                   break;
           }
       }

       DeletePixelHistos(verbosityLevel);

       if (verbosityLevel > 2)
       {
           cout << endl << "...End of Pixel"
                << endl << "------------------------------------------------"
                << endl;
       }
    }
    // End of Loop over all Pixels

//-------------------------------------
// Draw Histograms
//-------------------------------------

    SaveHistograms(     //save histograms of all pixel into output root file
                OutRootFileName,
                CreateSubDirName("All"),
                hAllPixelList,
                verbosityLevel
                );

//    SaveHistograms(     //save histograms of generell results into output root file
//                OutRootFileName,
//                "root",
//                hRootList,
//                verbosityLevel
//                );



    if (ProduceGraphic)
    {
        UpdateCanvases(
                    verbosityLevel,
                    MAX_PULS_ORDER
                    );
    }

    WriteAllPixelTemplateToCsv(
                OutPutPath,
                "PulseTemplate_PointSet",
                "Maximum",
                verbosityLevel
                );

    DeletePixelCanvases( verbosityLevel );
    return( 0 );
}
//----------------------------------------------------------------------------
// end of main function
//-----------------------------------------------------------------------------




//-----------------------------------------------------------------------------
// Funktions
//-----------------------------------------------------------------------------


void
BookPixelHistos( int verbosityLevel)
{
    if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
    hList = new TObjArray;
    TString histoname;
    for (int order = 0; order < MAX_PULS_ORDER; order ++)
    {
//        SetHistogramAttributes(
//                    "PeakOverlay",
//                    order,
//                    gOverlayAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hPixelOverlay";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hPixelOverlay[order] = new TH2F(
                    histoname,
                    "Overlay of detected pulses of one pulse order for one Pixel",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5 ,
                    512,
                    -55.5,
                    200.5
                    );
        hPixelOverlay[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        //hPixelProfile->SetBit(TH2F::kCanRebin);
        hList->Add( hPixelOverlay[order] );

    //------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "PeakMaximum",
//                    order,
//                    gMaximumAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hPixelMax";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hPixelMax[order] = new TH1F(
                    histoname,
                    "Maximum value of each slice in overlay plot",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5
                    );
        hPixelMax[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        hList->Add( hPixelMax[order] );

//------------------------------------------------------------------------
    //        SetHistogramAttributes(
    //                    "PeakMaximum",
    //                    order,
    //                    gMaximumAttrib,
    //                    MaxAmplOfFirstPulse
    //                    );
            histoname = "hPixelMedian";
            histoname += order;
            if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
            hPixelMedian[order] = new TH1F(
                        histoname,
                        "Median of each slice in overlay plot",
                        gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                        (-1*gPixelOverlayXaxisLeft)-0.5,
                        gPixelOverlayXaxisRight-0.5
                        );
            hPixelMedian[order]->SetAxisRange(
                        gBSLMean - 5,
                        (gGainMean*(order+1)) + 10,
                        "Y");
            hPixelMedian[order]->SetLineColor(kRed);
            hPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
            hPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
            hList->Add( hPixelMedian[order] );

//------------------------------------------------------------------------
    //        SetHistogramAttributes(
    //                    "PeakMaximum",
    //                    order,
    //                    gMaximumAttrib,
    //                    MaxAmplOfFirstPulse
    //                    );
            histoname = "hPixelMean";
            histoname += order;
            if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
            hPixelMean[order] = new TH1F(
                        histoname,
                        "Mean of each slice in overlay plot",
                        gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                        (-1*gPixelOverlayXaxisLeft)-0.5,
                        gPixelOverlayXaxisRight-0.5
                        );
            hPixelMean[order]->SetAxisRange(
                        gBSLMean - 5,
                        (gGainMean*(order+1)) + 10,
                        "Y");
            hPixelMean[order]->SetLineColor(kBlue);
            hPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
            hPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
            hList->Add( hPixelMean[order] );

//------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "PeakMaxGaus",
//                    order,
//                    gMaxGausAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hPixelEdgeOverlay";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hPixelEdgeOverlay[order] = new TH2F(
                    histoname,
                    "Overlay at rising edge of detected pulses of one pulse order for one Pixel ",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5 ,
                    512,
                    -55.5,
                    200.5
                    );

        hPixelEdgeOverlay[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        hList->Add( hPixelEdgeOverlay[order] );

    //------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "PeakProfile",
//                    order,
//                    gProfileAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hPixelProfile";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hPixelProfile[order] = new TProfile(
                    histoname,
                    "Mean value of each slice in overlay plot (Tprofile)",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
                    (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
                    gPixelOverlayXaxisRight-0.5 ,                    //xup
    //                512,                                            //nbinsy
                    -55.5,                                          //ylow
                    300.5,                                          //yup
                    "s");                                           //option
        hPixelProfile[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        //hPixelProfile->SetBit(TH2F::kCanRebin);
        hList->Add( hPixelProfile[order] );


    //------------------------------------------------------------------------
    //        SetHistogramAttributes(
    //                    "PeakProfile",
    //                    order,
    //                    gProfileAttrib,
    //                    MaxAmplOfFirstPulse
    //                    );
            histoname = "hPixelProfile2";
            histoname += order;
            if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
            hPixelProfile2[order] = new TProfile(
                        histoname,
                        "Mean value of each slice in overlay plot (Tprofile)",
                        gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
                        (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
                        gPixelOverlayXaxisRight-0.5 ,                    //xup
        //                512,                                            //nbinsy
                        -55.5,                                          //ylow
                        300.5,                                          //yup
                        "s");                                           //option
            hPixelProfile2[order]->SetLineColor(kRed);
            hPixelProfile2[order]->SetAxisRange(
                        gBSLMean - 5,
                        (gGainMean*(order+1)) + 10,
                        "Y");
            hPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
            hPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
            //hPixelProfile->SetBit(TH2F::kCanRebin);




    }
    if (verbosityLevel > 2) cout << "...done" << endl;
}
//end of BookPixelHistos
//----------------------------------------------------------------------------


void DeletePixelHistos( int verbosityLevel )
{
    if (verbosityLevel > 2) cout << endl << "...delete current pixel histograms" ;
    for (int order = 0; order < MAX_PULS_ORDER; order ++)
    {
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelOverlay" << order;
        delete hPixelOverlay[order];
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order;
        delete hPixelMax[order];
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelMedian" << order;
        delete hPixelMedian[order];
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelMean" << order;
        delete hPixelMean[order];
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order;
        delete hPixelEdgeOverlay[order];
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order;
        delete hPixelProfile[order];
        if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order;
        delete hPixelProfile2[order];
    }
    if (verbosityLevel > 3) cout << endl << "...deleting hList";
    delete hList;
    if (verbosityLevel > 2) cout << endl << "...done" << endl;
}
//end of DeletePixelHistos
//----------------------------------------------------------------------------

void BookTestHistos( int verbosityLevel )
{
    if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;

    hTesthisto = new TH1F (
                "hTesthisto",
                "Deviation of rising edge and maximum",
                600,
                -10.1,
                10.1
                );
    hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    hTesthisto->GetYaxis()->SetTitle( "counts" );

    hTesthisto2 = new TH2F (
                "hTesthisto2",
                "Deviation of rising edge and maximum by event #",
                gNEvents,
                250,
                800,
                600,
                -10.1,
                10.1
                );
//    hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
//    hTesthisto2->GetYaxis()->SetTitle( "counts" );
//    hTesthisto2->SetMarkerStyle( 2 );

    if (verbosityLevel > 2) cout << "...done" << endl;
}
//end of BookTestHistos
//----------------------------------------------------------------------------

void BookHistos( int verbosityLevel )
{
        if (verbosityLevel > 2) cout << endl << "...book histograms" << endl;
    hAllPixelList = new TObjArray;
    hRootList = new TObjArray;
    debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ];
    for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){
        debugHistos[ type ].SetBins(1024, 0, 1024);
        debugHistos[ type ].SetLineColor(1);
        debugHistos[ type ].SetLineWidth(2);
        debugHistos[ type ].SetStats(false);

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

        // set Y axis paras
        debugHistos[ type ].GetYaxis()->SetLabelSize(0.08);
        debugHistos[ type ].GetYaxis()->SetTitleSize(0.08);
        debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
        debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]");
    }

    // All Pixel Histograms
    //------------------------------------------------------------------------
    TString histoname;
    for (int order = 0; order < MAX_PULS_ORDER; order ++)
    {
//        SetHistogramAttributes(
//                    "AllPixelPeakOverlay",
//                    order,
//                    gOverlayAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelOverlay";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelOverlay[order] = new TH2F(
                    histoname,
                    "Overlay of detected Pulses of one Order for All Pixels",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5 ,
                    512,
                    -55.5,
                    300.5
                    );
        hAllPixelOverlay[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        //hAllPixelOverlay->SetBit(TH2F::kCanRebin);
        hAllPixelList->Add( hAllPixelOverlay[order] );

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

//        SetHistogramAttributes(
//                    "AllPixelPeakMax",
//                    order,
//                    gMaximumAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelMax";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelMax[order] = new TH1F(
                    histoname,
                    "Maximum value of each slice in overlay plot for All Pixels",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5
                    );
        hAllPixelMax[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        hAllPixelList->Add( hAllPixelMax[order] );

//------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "PeakMaximum",
//                    order,
//                    gMaximumAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelMedian";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelMedian[order] = new TH1F(
                    histoname,
                    "Median of each slice in overlay plot for All Pixels",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5
                    );
        hAllPixelMedian[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelMedian[order]->SetLineColor(kRed);
        hAllPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        hAllPixelList->Add( hAllPixelMedian[order] );
        if (verbosityLevel > 3) cout << "...BREAKPOINT in " << histoname << endl;
//------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "PeakMaximum",
//                    order,
//                    gMaximumAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelMean";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelMean[order] = new TH1F(
                    histoname,
                    "Mean of each slice in overlay plot for All Pixels",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5
                    );
        hAllPixelMean[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelMean[order]->SetLineColor(kBlue);
        hAllPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        hAllPixelList->Add( hAllPixelMean[order] );

//------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "PeakMaxGaus",
//                    order,
//                    gMaxGausAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelEdgeOverlay";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelEdgeOverlay[order] = new TH2F(
                    histoname,
                    "Overlay at rising edge of detected pulses of one pulse order for All Pixels ",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
                    (-1*gPixelOverlayXaxisLeft)-0.5,
                    gPixelOverlayXaxisRight-0.5 ,
                    512,
                    -55.5,
                    200.5
                    );

        hAllPixelEdgeOverlay[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        hAllPixelList->Add( hAllPixelEdgeOverlay[order] );

//------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "AllPixelPeakProfile",
//                    order,
//                    gProfileAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelProfile";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelProfile[order] = new TProfile(
                    histoname,
                    "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
                    (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
                    gPixelOverlayXaxisRight-0.5 ,                    //xup
    //                512,                                            //nbinsy
                    -55.5,                                          //ylow
                    300.5,                                          //yup
                    "s");                                           //option
        hAllPixelProfile[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        //hPixelProfile->SetBit(TH2F::kCanRebin);
        hAllPixelList->Add( hAllPixelProfile[order] );

//------------------------------------------------------------------------
//        SetHistogramAttributes(
//                    "AllPixelPeakProfile",
//                    order,
//                    gProfileAttrib,
//                    MaxAmplOfFirstPulse
//                    );
        histoname = "hAllPixelProfile2";
        histoname += order;
        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
        hAllPixelProfile2[order] = new TProfile(
                    histoname,
                    "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
                    (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
                    gPixelOverlayXaxisRight-0.5 ,                    //xup
    //                512,                                            //nbinsy
                    -55.5,                                          //ylow
                    300.5,                                          //yup
                    "s");                                           //option
        hAllPixelProfile2[order]->SetAxisRange(
                    gBSLMean - 5,
                    (gGainMean*(order+1)) + 10,
                    "Y");
        hAllPixelProfile2[order]->SetLineColor(kRed);
        hAllPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
        hAllPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
        //hPixelProfile->SetBit(TH2F::kCanRebin);
        hAllPixelList->Add( hAllPixelProfile2[order] );



    }//end of for over order
    if (verbosityLevel > 2) cout << "...done" << endl;
}
// end of BookHistos
//----------------------------------------------------------------------------


void CreateRootFile(const char * loc_fname, int verbosityLevel)
{
    TFile* tf = new TFile(loc_fname,"UPDATE");
    if (tf->IsZombie())
    {
        cout << "Error opening file" << endl;
        exit(-1);
    }
    else {
        if (verbosityLevel > 1) cout << "creating root-file successfull" << endl;
        tf->Close();
    }
}
//end of CreateRootFile
//----------------------------------------------------------------------------


TFile* OpenRootFile(const char * loc_fname, int verbosityLevel)
{
    TFile* tf = new TFile(loc_fname,"UPDATE");
    if (tf->IsZombie())
    {
        cout << "Error opening file" << endl;
        exit(-1);
    }
    else {
        if (verbosityLevel > 1) cout << "...opening root-file:"
                                     << loc_fname
                                     << " successfull" << endl;
    }
    return tf;
}
//end of OpenRootFile
//----------------------------------------------------------------------------


void CloseRootFile(TFile* tf)
{
    if (tf->IsZombie())
    {
    cout << "Error closing file" << endl;
    exit(-1);
    } else {
    tf->Close();
    }
}
//end of CloseRootFile
//----------------------------------------------------------------------------


void SaveHistograms(
        const char* loc_fname,
        const char* subdirectory,
        TObjArray* histList,
        int verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...saving pixel histograms to subdirectory: "
                                 << subdirectory << endl ;

    TFile* tf = OpenRootFile(loc_fname, verbosityLevel);
    if ( subdirectory != "root")
    {
        tf->cd();
        tf->mkdir(subdirectory,subdirectory);
        tf->cd(subdirectory);
    }
    else
    {
        tf->cd();
    }
    histList->Write(); // write the major histograms into the top level directory
    if (verbosityLevel > 3) tf->ls();
    CloseRootFile( tf ); // close the file
    if (verbosityLevel > 2) cout << "...done" << endl;
}
// end of SaveHistograms
//----------------------------------------------------------------------------


TString CreateSubDirName(int pixel)
{
    TString path = "Pixel_";
    path += pixel;
//    cout << "path " << path << endl ;
    return path;
}
// end of CreateSubDirName
//----------------------------------------------------------------------------

TString CreateSubDirName(const char* title)
{
    TString path = title;
    path += "_Pixel";
//    cout << "path " << path << endl ;
    return path;
}
// end of CreateSubDirName
//----------------------------------------------------------------------------

void FillHistograms(
        vector<Region>* pZXings,
        int AmplWindowWidth,
        int eventNumber,
        int verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ;
    bool use_this_peak=false;
    int order_of_pulse=0;
    vector<Region>::iterator reg;

    //Loop over all found zerocrossings in Region
    for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
    {
        //skip those who are at the Rim of the ROI
        if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12)
        {
            if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl;
            continue;
        }
        //domstest->Fill(reg->maxPos);

        // Define axis range of Histogram
        int Left = reg->maxPos - gPixelOverlayXaxisLeft;
        int Right = reg->maxPos + gPixelOverlayXaxisRight;
        int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft;
        int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight;

        if (Left < 0) Left =0;
        if (Right > (int)Ameas.size() ) Right=Ameas.size() ;
        if (EdgeLeft < 0) EdgeLeft =0;
        if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ;


        hTesthisto->Fill( reg->slopeOfRisingEdge ) ;
        hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ;

        if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;

        //determine order of pulse and which histogram shall be filled
        if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
        for(int order = 0; order < MAX_PULS_ORDER; order++)
        {
            if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2))
                && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2))
            {
                //hOverlayTemp = hPixelOverlay[order];
                if (verbosityLevel > 2) cout << "...#" << order ;
                order_of_pulse = order;
                use_this_peak = true;
                break;
            }
         else if (order >= (MAX_PULS_ORDER - 1)){
                use_this_peak = false;
                if (verbosityLevel > 2) cout << "...NONE" << endl ;
         }
        }

        //Fill overlay und profile histograms
        if (use_this_peak){
            if (verbosityLevel > 2) cout << "...filling" ;
            for ( int pos = Left; pos < Right; pos++){
//                if ();
                hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
                hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
                hAllPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
                hAllPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
            }
            for ( int pos = EdgeLeft; pos < EdgeRight; pos++){
//                cout << endl << "###filling edge histo###" << endl;
                hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
                hPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
                hAllPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
                hAllPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
//                cout << endl << "######" << endl;
            }
            if (verbosityLevel > 2) cout << "...done" << endl;
        }
        //Breakout option
        if (breakout)
        break;
   }
    if (verbosityLevel > 2) cout << "...done" << endl;
}
// end of FillHistograms
//----------------------------------------------------------------------------

void DrawTestHistograms(
        int verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ;

    cgpTestHistos->cd(1);
    hTesthisto->Draw();
    cgpTestHistos->cd(2);
    hTesthisto2->Draw("BOX");
}
// end of DrawTestHistograms
//----------------------------------------------------------------------------

void DrawPulseHistograms(
        int verbosityLevel,
        int max_pulse_order
        )
{
    if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ;
    for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
    {
        cgpPixelPulses[pulse_order]->cd(1);
        hPixelOverlay[pulse_order]->Draw("COLZ");
        cgpPixelPulses[pulse_order]->cd(3);
        hPixelProfile[pulse_order]->Draw("COLZ");
        hPixelProfile2[pulse_order]->Draw("SAME");

        cgpPixelPulses[pulse_order]->cd(2);
        hPixelEdgeOverlay[pulse_order]->Draw("COLZ");
//        cgpPixelPulses[pulse_order]->cd(4);
//        hTesthisto->Draw();

        cgpAllPixelPulses[pulse_order]->cd(1);
        hAllPixelOverlay[pulse_order]->Draw("COLZ");
        cgpAllPixelPulses[pulse_order]->cd(2);
        hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ");

        cgpAllPixelPulses[pulse_order]->cd(3);
        hAllPixelProfile[pulse_order]->Draw("COLZ");
        hAllPixelProfile2[pulse_order]->Draw("SAME");

    }
}
// end of DrawPulseHistograms
//----------------------------------------------------------------------------


void DrawMaximumHistograms(
        int verbosityLevel,
        int max_pulse_order
        )
{
    if (verbosityLevel > 2) cout << endl << "...drawing maximum histograms" ;
    for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
    {
        cgpPixelPulses[pulse_order]->cd(4);
        hPixelMax[pulse_order]->Draw();
        hPixelMedian[pulse_order]->Draw("SAME");
        hPixelMean[pulse_order]->Draw("SAME");

        cgpAllPixelPulses[pulse_order]->cd(4);
        hAllPixelMax[pulse_order]->Draw();
        hAllPixelMedian[pulse_order]->Draw("SAME");
        hAllPixelMean[pulse_order]->Draw("SAME");
//        cgpAllPixelPulses[pulse_order]->cd(3);
     //   hAllPixelMaxGaus[pulse_order]->Draw("COLZ");
    }
}
// end of DrawMaximumHistograms
//----------------------------------------------------------------------------


void UpdateCanvases(
        int verbosityLevel,
        int max_pulse_order
        )
{
    if (verbosityLevel > 3) cout << endl << "...updating canvases" ;
    for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
    {
        cgpPixelPulses[pulse_order]->Modified();
        cgpPixelPulses[pulse_order]->Update();
        cgpAllPixelPulses[pulse_order]->Modified();
        cgpAllPixelPulses[pulse_order]->Update();
        cgpTestHistos->Modified();
        cgpTestHistos->Update();
    }
}
// end of UpdateCanvases
//----------------------------------------------------------------------------



void
PlotMedianEachSliceOfPulse(
        unsigned int    max_pulse_order,
        int             fitdata,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of slice's Median"
                                 << endl;
    for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    {
        if (verbosityLevel > 2) cout << "\t...calculation of "
                                     << "pulse order " << pulse_order;

        Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();

        for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
           TH1 *h1 = hPixelOverlay[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
           float median = MedianOfH1(h1);
           float mean = h1->GetMean();
           if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
           delete h1;
           hPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
           hPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
           hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
           hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
        }

        if (verbosityLevel > 2) cout << "\t...done" << endl;

        if (fitdata)
        {
           FitMaxPropabilityPuls(
                    hPixelMedian[pulse_order],
                    verbosityLevel
                    );
        }
    }
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------

Double_t MedianOfH1 (
        TH1*            h1
        )
{
   //compute the median for 1-d histogram h1
   Int_t nbins = h1->GetXaxis()->GetNbins();
   Double_t *x = new Double_t[nbins];
   Double_t *y = new Double_t[nbins];
   for (Int_t i=0;i<nbins;i++) {
      x[i] = h1->GetXaxis()->GetBinCenter(i+1);
      y[i] = h1->GetBinContent(i+1);
   }
   Double_t median = TMath::Median(nbins,x,y);
   delete [] x;
   delete [] y;
   return median;
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------

void
PlotPulseShapeOfMaxPropability(
        unsigned int    max_pulse_order,
        int             fitdata,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of max. propability"
                                 << endl;
    for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    {
        if (verbosityLevel > 2) cout << "\t...calculation of "
                                     << "pulse order " << pulse_order;
        //  vector max_value_of to number of timeslices in Overlay Spectra
        float max_value_of_slice;

        Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
        if (verbosityLevel > 2) cout << "...generating projections" << endl;
        for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
        {
            if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;
            //put maximumvalue of every Y-projection of every timeslice into vector
            TH1D *h1 =	hPixelOverlay[pulse_order]->ProjectionY( "", TimeSlice,TimeSlice);

            max_value_of_slice = h1->GetBinCenter( h1->GetMaximumBin() );

            if (verbosityLevel > 2) cout << " with max value "
                                         << max_value_of_slice << endl;
            hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
            hAllPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
            delete h1;
        }

        if (verbosityLevel > 2) cout << "\t...done" << endl;

        if (fitdata)
        {
           FitMaxPropabilityPuls(
                    hPixelMax[pulse_order],
                    verbosityLevel
                    );
        }
    }
}


void
FitMaxPropabilityPuls(
        TH1F*       hMaximumTemp,
        int         verbosityLevel
        )
    {
      if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
      if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
      hMaximumTemp->Fit("landau", "", "", -50, 250);
      if (verbosityLevel > 2) cout << "...done" << endl;
    }


//void SetHistogrammSettings(
//        const char*     histoType,
//        int             max_puls_order,
//        int             MaxAmplOfFirstPulse,
//        int             verbosityLevel
//        )
//{
//    if (verbosityLevel > 2) cout << endl << "...set histogram Settings" << endl;
//    for (int pulse_order = 1; pulse_order <= max_puls_order; pulse_order++)
//    {
//        qHistSetting[pulse_order].yMax = (MaxAmplOfFirstPulse * pulse_order);
//        cout << "\t Max @ " << gHistSetting[pulse_order].yMax << endl;
//        gHistSetting[pulse_order].name = "h";
//        gHistSetting[pulse_order].name += histoType;
//        gHistSetting[pulse_order].name += "PeakOverlay";
//        gHistSetting[pulse_order].name += pulse_order;
//        gHistSetting[pulse_order].title = "PeakOverlay with Amplitude around ";
//        gHistSetting[pulse_order].title += gHistSetting[pulse_order].yMax;
//        gHistSetting[pulse_order].title += " mV";
//        gHistSetting[pulse_order].Mname = "h";
//        gHistSetting[pulse_order].Mname += histoType;
//        gHistSetting[pulse_order].Mname += "PeakMaximum";
//        gHistSetting[pulse_order].Mname += pulse_order;
//        gHistSetting[pulse_order].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of ";
//        gHistSetting[pulse_order].Mtitle += gHistSetting[pulse_order].yMax;
//        gHistSetting[pulse_order].Mtitle += " mV";
//    }
//    if (verbosityLevel > 2) cout << "...done" << endl;
//}


//void CreatePulseCanvas(
//        TCanvas*    array_of_canvases,
//        unsigned int        max_pulse_order,
//        const char*      canvas_name,
//        const char*      canvas_title,
//        const char*      pixel,
//        int verbosityLevel
//        )
//{
//    if (verbosityLevel > 2) cout << "...building Canvases" ;

//    for (
//         unsigned int pulse_order = 1;
//         pulse_order <= max_pulse_order;
//         pulse_order++
//          )
//    {
//        TString name = canvas_name;
//                name += " ";
//                name += pulse_order;
//                name += " ";

//                TString title = "Order ";
//                title += pulse_order;
//                title += " ";
//                title += canvas_title;
//                title += " of ";
//                title += pixel;
//                title += ". Pixel";

//        array_of_canvases[pulse_order] = new TCanvas(name, title, 1, 1, 1400, 700);
//        array_of_canvases[pulse_order].Divide(2,2);
//    }
//    if (verbosityLevel > 2) cout << "...done" << endl;
//}


//void
//DeletePulseCanvas(
//        TCanvas* array_of_canvases,
//        unsigned int     max_pulse_order
//        )
//{
//    for (
//         unsigned int pulse_order = 1;
//         pulse_order <= max_pulse_order;
//         pulse_order++
//          )
//    {
//        delete array_of_canvases[pulse_order];
//    }
//}




void
DeletePixelCanvases( int verbosityLevel )
{
    if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ;
    for (int pulse_order = 0; pulse_order < MAX_PULS_ORDER; pulse_order ++)
    {
        delete cgpPixelPulses[pulse_order];
    }
}

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

//void SetHistogramAttributes(
//        const char*     histoType,
//        int             pulse_order,
//        histAttrib_t*   histoAttrib,
//        int             max_ampl_of_first_pulse
//        )
//{
//    histoAttrib[pulse_order].yMax = (max_ampl_of_first_pulse * pulse_order);
//    histoAttrib[pulse_order].yMin = (histoAttrib[pulse_order].yMax/2);
//    histoAttrib[pulse_order].name = "h";
//    histoAttrib[pulse_order].name += histoType;
//    histoAttrib[pulse_order].name = "_";
//    histoAttrib[pulse_order].name += pulse_order;

//    histoAttrib[pulse_order].title += histoType;
//    histoAttrib[pulse_order].title += " with Amplitude around ";
//    histoAttrib[pulse_order].title += histoAttrib[pulse_order].yMax;
//    histoAttrib[pulse_order].title += " mV";

//    histoAttrib[pulse_order].xTitle += "Timeslices [a.u.]";
//    histoAttrib[pulse_order].yTitle += "Amplitude [mV]";
//}

void WritePixelTemplateToCsv(
        TString         path,
        const char*     csv_file_name,
        const char*     overlay_method,
        int             pixel,
        int             verbosityLevel
        )
{
    path += csv_file_name;
    path += "_";
    path += pixel;
    path += ".csv";

    Int_t nbins = hPixelMax[0]->GetXaxis()->GetNbins();
    if (verbosityLevel > 0)
    {
        cout << "writing point-set to csv file: " ;
        cout << path << endl;
        cout << "...opening file" << endl;
    }
    if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
    ofstream out;
    out.open( path );
    out << "### point-set of a single photon pulse template" << endl;
    out << "### template determined with pulse overlay at: "
        << overlay_method << endl;
    out << "### Slice's Amplitude determined by calculating the " << endl
        << "### value of maximum propability of slice -> Amplitude1 " << endl
        << "### mean of slice -> Amplitude2 " << endl
        << "### median of slice -> Amplitude3 " << endl
        << "### for each slice" << endl;
    out << "### Pixel number (CHid): " << pixel << endl
        << endl;

    out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;

    for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    {
        out << TimeSlice << "," ;
        out << hPixelMax[0]->GetBinContent(TimeSlice) << ",";
        out << hPixelMean[0]->GetBinContent(TimeSlice) << ",";
        out << hPixelMedian[0]->GetBinContent(TimeSlice) << endl;
    }
    out.close();
    if (verbosityLevel > 0) cout << "...file closed" << endl;
}

void WriteAllPixelTemplateToCsv(
        TString         path,
        const char*     csv_file_name,
        const char*     overlay_method,
        int             verbosityLevel
        )
{
    path += csv_file_name;
    path += "_AllPixel.csv";

    Int_t nbins = hAllPixelMax[0]->GetXaxis()->GetNbins();
    if (verbosityLevel > 0)
    {
        cout << "writing point-set to csv file: " ;
        cout << path << endl;
        cout << "...opening file" << endl;
    }
    if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
    ofstream out;
    out.open( path );
    out << "### point-set of a single photon pulse template" << endl;
    out << "### template determined with pulse overlay at: "
        << overlay_method << "of all Pixels";
    out << "### Slice's Amplitude determined by calculating the " << endl
        << "### value of maximum propability of slice -> Amplitude1 " << endl
        << "### mean of slice -> Amplitude2 " << endl
        << "### median of slice -> Amplitude3 " << endl
        << "### for each slice" << endl
        << endl;

    out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;

    for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    {
        out << TimeSlice << "," ;
        out << hAllPixelMax[0]->GetBinContent(TimeSlice) << ",";
        out << hAllPixelMean[0]->GetBinContent(TimeSlice) << ",";
        out << hAllPixelMedian[0]->GetBinContent(TimeSlice) << endl;
    }
    out.close();
    if (verbosityLevel > 0) cout << "...file closed" << endl;
}