/********************** 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 FCalcPulseTemplate(
  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           sampleSetSize       = 36,
  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 Pixel Sets
//-------------------------------------
    for ( int firstPixelOfSample = firstpixel;
          firstPixelOfSample < firstpixel + npixel;
          firstPixelOfSample += sampleSetSize )
    {

        if (verbosityLevel >= 0)
        {
            cout << "------------------------------------------------" << endl
                 << "...processing Sample from Pixel: "
                 << firstPixelOfSample;
                 << "to Pixel: "
                 << firstPixelOfSample+sampleSetSize-1 << 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;
            }

//--------------------------------------------------------------------
// Loops over every Pixel of a Set of Pixels
//--------------------------------------------------------------------
        for ( int pixel = firstPixelOfSample;
            pixel < sampleSetSize || pixel < firstpixel + npixel;
            pixel++ )
        {
            if (verbosityLevel > 0)
           {
             cout << "-------------------------------------" << endl
                  << "...processing Pixel: " << pixel
                  << "/" << nPixelSample << 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 Sample
//-------------------------------------
// end of Loops over Sample
//-------------------------------------

        }//End of Loop over Events
//-------------------------------------
// end of Loops over Events
//-------------------------------------




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

    PlotPulseShapeOfMaxPropability(
            MAX_PULS_ORDER,
            fitdata,
            verbosityLevel
         );

    for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    {
    if (verbosityLevel > 2) cout << "\t...plotting most probable pulse of"
                                 << "pulse order " << pulse_order;
    PlotMedianEachSliceOfPulse(
         hPixelOverlay[pulse_order],
         hPixelMedian[pulse_order],
         fitdata,
         verbosityLevel
         );

    PlotMeanEachSliceOfPulse(
         hPixelOverlay[pulse_order],
         hPixelMean[pulse_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
//-----------------------------------------------------------------------------

/*
ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel
wie viel speicher braucht das?
Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren
am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute.
ne klasse wäre super

*/
void
BookPixelHistos(
        int     sampleSize,

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