#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <iostream>

#include <TCanvas.h>
#include <TTimer.h>
#include <Getline.h>

#include "templateextractors.h"

using namespace std;
using namespace TMath;


void
CalcMaxPropabilityOfSlice(
        TH2*            inputHisto,
        TH1*            outputHisto,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2)
        cout << endl
             << "...calculating slieces value of max. propability"
             << endl;

    float   sliceMax    = 0;
    TH1D*   hTemp       = NULL;
    int     nbins       = inputHisto->GetXaxis()->GetNbins();

    if (verbosityLevel > 2)
        cout << "...generating projections" << endl;

    for (Int_t slice = 1; slice <= nbins; slice++)
    {
        if (verbosityLevel > 2)
            cout << "...slice: " << slice;

    //put maximumvalue of every Y-projection of every slice into vector
        hTemp       = inputHisto->ProjectionY( "", slice,slice);
        sliceMax    = hTemp->GetBinCenter( hTemp->GetMaximumBin() );

        if (verbosityLevel > 2)
            cout << " with max value "
                 << sliceMax << endl;

        outputHisto->SetBinContent(slice, sliceMax );
    }

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


void
PlotMaxPropabilityPulse(
        Pixel*          CurrentPixel,
        TString         overlayType,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of max. propability"
                                 << endl;
    TH2F**  hInputHistoArray    = NULL;
    TH1F**  hOutputHistoArray   = NULL;

    if (overlayType.Contains("Edge"))
    {
        hInputHistoArray    = CurrentPixel->hEdgeOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelEdgeMax;

    }
    else if (overlayType.Contains("Maximum"))
    {
        hInputHistoArray    = CurrentPixel->hMaxOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelMax;
    }
    else
    {
        cout << endl << "Unknown Overlay Method-->aborting" << endl;
        return;
    }

    for ( int pulse_order = 0 ;
          pulse_order < CurrentPixel->mMaxPulseOrder ;
          pulse_order ++)
    {
        if (verbosityLevel > 2) cout << "\t...calculation of "
                                     << "pulse order " << pulse_order;
        //  vector max_value_of to number of slices in Overlay Spectra
        CalcMaxPropabilityOfSlice(
                    hInputHistoArray[pulse_order],
                    hOutputHistoArray[pulse_order],
                    verbosityLevel);

        CalcMaxPropabilityOfSlice(
                    hInputHistoArray[pulse_order],
                    hOutputHistoArray[pulse_order],
                    verbosityLevel);
    }
}
// end of PlotMaxPropabilityPulse
//----------------------------------------------------------------------------

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

double MedianOfH1 ( //THE PROBLEM MUST BE HERE!!!!!!!
        TH1D*            inputHisto
        )
{
   //compute the median for 1-d histogram h1
   Int_t nbins = inputHisto->GetXaxis()->GetNbins();
   double *x = new double[nbins];
   double *y = new double[nbins];
   for (Int_t i=0;i<nbins;i++) {
      x[i] = inputHisto->GetXaxis()->GetBinCenter(i);
      y[i] = inputHisto->GetBinContent(i);
   }
   double median = TMath::Median(nbins,x,y);
   delete [] x;
   delete [] y;
   return median;
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------

Double_t MedianOfH1withRndSlices (
        TH1*            inputHisto,     //histogram from which median will be calculated
        vector<int>*    position        //array with random slice numbers
        )
{
   //compute the median for 1-d histogram h1
   Int_t nbins = inputHisto->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] = inputHisto->GetXaxis()->GetBinCenter(position->at(i) );
      y[i] = inputHisto->GetBinContent(position->at(i) );
   }
   Double_t median = TMath::Median(nbins,x,y);
   delete [] x;
   delete [] y;
   return median;
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------

int ExtractTH1EnriesToVector (
        TH1*            inputHisto,     //histogram from which median will be calculated
        vector<int>*    entries        //array with random slice numbers
        )
{
   //compute number of bins in imput histo
   Int_t nbins = inputHisto->GetXaxis()->GetNbins();
   int quantity;
   for (Int_t i=0;i<nbins;i++) {
       quantity = inputHisto->GetBinContent(i);
       for (int j = 0; j < quantity; j++){
           entries->push_back(inputHisto->GetXaxis()->GetBinCenter(i));
       }
   }
   return entries->size();
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------

void CalculateErrorsWithBootstrapping(TH1* inputHisto, int numIt, double* par, double* parErr)
{
//    cout << "...calculating errors of " << inputHisto->GetName() << endl;
    Sample  sample;
    double  Mean[numIt];
    double  Median[numIt];
    double  Max[numIt];
    double  parameter[3]    = {0,0,0};
    double  parameterErr[3] = {0,0,0};
    for (int i = 0; i < numIt; i++)
    {
        Mean[i]     = 0;
        Median[i]   = 0;
        Max[i]      = 0;
    }

//    TCanvas test_canvas;
//    test_canvas.Divide(2,1);
//    test_canvas.cd(1);
//    inputHisto->Draw();

    for (int i = 0; i < numIt; i++)
    {
//        test_canvas.cd(2);
        TH1D* tempHisto = (TH1D*)inputHisto->Clone("tempHisto");
//                    "",
//                    "",
//                    inputHisto->GetNbinsX(),
//                    inputHisto->GetXaxis()->GetBinLowEdge(  inputHisto->GetXaxis()->GetFirst() ),
//                    inputHisto->GetXaxis()->GetBinUpEdge(   inputHisto->GetXaxis()->GetLast() )
//                    );
        tempHisto->Reset();
        sample.SetSeed(sample.GetSeed()+1);
        sample.BootstrapTH1(inputHisto, tempHisto);

        Mean[i]     = tempHisto->GetMean();
        Median[i]   = MedianOfH1 ( tempHisto );
        Max[i]      = tempHisto->GetBinCenter( tempHisto->GetMaximumBin() );

        //improved determination of modus
//        TF1 gaus("fgaus", "gaus", Max[i]-10, Max[i]+10);
//        tempHisto->Fit("fgaus", "WWQRN0");
//        Max[i]      = gaus.GetParameter(1);


//        sample.SetSeed(sample.GetSeed() + 1);

        /*
        //Process gui events asynchronously during input
        tempHisto->Draw();
        test_canvas.Modified();
        test_canvas.Update();
        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")
        {
            return;
        }
        */



        delete tempHisto;
    }

    parameter[0]    = TMath::Mean(  numIt, Max);
    parameterErr[0] = TMath::RMS(   numIt, Max);
    parameter[1]    = TMath::Mean(  numIt, Median);
    parameterErr[1] = TMath::RMS(   numIt, Median);
    parameter[2]    = TMath::Mean(  numIt, Mean);
    parameterErr[2] = TMath::RMS(   numIt, Mean);

    for (int i = 0; i < 3; i++)
    {
        if (&par[i] != NULL)
        {
            par[i] =  parameter[i];
        }
        else cout << "par[" << i << "] array to small" << endl;
        if (&par[i] != NULL)
        {
            parErr[i] =  parameterErr[i];
        }
        else cout << "parErr[" << i << "] array to small" << endl;
    }

    return;
}

void
PlotMedianOfSlice(
        Pixel*          CurrentPixel,
        TString         overlayType,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of slice's Median"
                                 << endl;

    TH2F**  hInputHistoArray    = NULL;
    TH1F**  hOutputHistoArray   = NULL;
//    TH1*    hTempHisto          = NULL;
//    float   median              = 0;

    if (overlayType.Contains("Edge"))
    {
        hInputHistoArray    = CurrentPixel->hEdgeOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelEdgeMean;

    }
    else if (overlayType.Contains("Maximum"))
    {
        hInputHistoArray    = CurrentPixel->hMaxOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelMedian;
    }
    else
    {
        cout << endl << "Unknown Overlay Method-->aborting" << endl;
        return;
    }

    for (int pulse_order = 0 ;
         pulse_order < CurrentPixel->mMaxPulseOrder ;
         pulse_order ++)
    {
        if (verbosityLevel > 2) cout << "\t...calculation of "
                                     << "pulse order " << pulse_order;

        MedianOfTH2Slices(
                hInputHistoArray[pulse_order],
                hOutputHistoArray[pulse_order],
                verbosityLevel
                );

        ErrorMedianOfTH2Slices(
                hInputHistoArray[pulse_order],
                hOutputHistoArray[pulse_order],
                5,  //numIterations
                verbosityLevel
                );

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

//        for (Int_t slice=1;slice<=nbins;slice++) {

//            hTempHisto  = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
//            median      = MedianOfH1(hTempHisto);

//            if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);

//            hOutputHistoArray[pulse_order]->SetBinContent(slice, median );
////            delete h1;
//        }

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

void
MedianOfTH2Slices(
        TH2*        inputHisto,
        TH1*        outputHisto,
        int         verbosityLevel
        )
{
    Int_t   nbins       = inputHisto->GetXaxis()->GetNbins();
    float   median              = 0;

    for (Int_t slice=1;slice<=nbins;slice++) {

        median      = MedianOfH1( inputHisto->ProjectionY("",slice,slice) );

        if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);

        outputHisto->SetBinContent(slice, median );
//            delete h1;
    }
    return;
}
// end of MedianOfTH2Slices
//----------------------------------------------------------------------------

double
ErrorMedianOfH1(
        TH1*        inputHisto,
        int         numIterations,
        int         verbosityLevel
        )
{
    if (verbosityLevel > 5) cout << "...calculating errors of median" << endl;
//    float MedianOfSliceMean;
    double medianRMS = 0;

    int sample_min  = inputHisto->GetXaxis()->GetFirst();
    int sample_max  = inputHisto->GetXaxis()->GetLast();
    int sample_size = inputHisto->GetXaxis()->GetNbins();
//    int sample_size = sample_max - sample_min;
//    cout << "sample_min " << sample_min
//         << ", sample_max " << sample_max
//         << ", sample_size " << sample_size
//         << endl;
    double   median[numIterations];
    vector<int> rndList;
    Sample sample(sample_size);

    for (int i = 0; i < numIterations; i++)
    {
        sample.SetSeed(sample.GetSeed() + 1);
        sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
        median[i] = MedianOfH1withRndSlices(
                    inputHisto,
                    &rndList
                    );
    }

//    MedianOfSliceMean    = TMath::Mean(numIterations, median);
    medianRMS     = RMS(numIterations, median);
//    if (verbosityLevel > 2) printf("Mean Median=%g\n",MedianOfSliceMean);
    if (verbosityLevel > 5) printf("medianRMS=%g\n",medianRMS);

    return medianRMS;
}

void
ErrorMedianOfTH2Slices(
        TH2*        inputHisto,
        TH1*        outputHisto,
        int         numIterations,
        int         verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << "...calculating errors of median" << endl;
    Int_t   nbins       = inputHisto->GetXaxis()->GetNbins();
    cout << "nbins " << nbins << endl;
    float MedianOfSliceMean[nbins];
    float MedianOfSliceRMS[nbins];
    float median[numIterations];
    vector<int> rndList;
//    rndList->clear();

    TH1 *htemp = NULL;

    for (Int_t slice=1;slice<=nbins;slice++) {
        htemp   = inputHisto->ProjectionY("",slice,slice);


        int sample_min  = htemp->GetXaxis()->GetFirst();
        int sample_max  = htemp->GetXaxis()->GetLast();
        int sample_size = htemp->GetXaxis()->GetNbins();

//        cout << "sample_min " << sample_min
//             << ", sample_max " << sample_max
//             << ", sample_size " << sample_size
//             << endl;

        Sample sample(sample_size);
        for (int i = 0; i < numIterations; i++)
        {
            sample.SetSeed(sample.GetSeed() + 1);
            sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
            median[i] = MedianOfH1withRndSlices(
                        htemp,
                        &rndList
                        );

        }
        MedianOfSliceMean[slice]    = TMath::Mean(numIterations, median);
        MedianOfSliceRMS[slice]     = TMath::RMS(numIterations, median);
        outputHisto->SetBinError(slice, MedianOfSliceRMS[slice]);
//        outputHisto->SetBinError(slice, RMS(numIterations, median) );
        if (verbosityLevel > 2) printf("Mean Median of Slice %d, Mean Median=%g\n",slice,MedianOfSliceMean[slice]);
        if (verbosityLevel > 2) printf("RMS of Median of Slice %d, MedianRMS=%g\n",slice,MedianOfSliceRMS[slice]);
//        outputHisto->SetBinContent(slice, median );
//            delete h1;
    }
    return;
}
// end of ErrorMedianOfTH2Slices
//----------------------------------------------------------------------------

void
PlotMeanOfSlice(
        Pixel*          CurrentPixel,
        TString         overlayType,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of slice's Mean"
                                 << endl;

    TH2F**  hInputHistoArray    = NULL;
    TH1F**  hOutputHistoArray   = NULL;
    TH1*    hTempHisto          = NULL;
    float   mean              = 0;

    if (overlayType.Contains("Edge"))
    {
        hInputHistoArray    = CurrentPixel->hEdgeOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelEdgeMean;

    }
    else if (overlayType.Contains("Maximum"))
    {
        hInputHistoArray    = CurrentPixel->hMaxOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelMean;
    }
    else
    {
        cout << endl << "Unknown Overlay Method-->aborting" << endl;
        return;
    }

    for (int pulse_order = 0 ;
         pulse_order < CurrentPixel->mMaxPulseOrder ;
         pulse_order ++)
    {
        if (verbosityLevel > 2) cout << "\t...calculation of "
                                     << "pulse order " << pulse_order;

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

        for (Int_t slice=1;slice<=nbins;slice++) {

            hTempHisto  = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
            mean      = hTempHisto->GetMean();

            if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",slice,mean);

            hOutputHistoArray[pulse_order]->SetBinContent(slice, mean );
//            delete h1;
        }

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

void
ExtractPulseTemplate(
        Pixel*          CurrentPixel,
        TString         overlayType,
        int             pulse_order,
        int             numBootstrapIt,// = 10, //number of iterations for bootstrapping
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape"
                                 << " of max. propability"
                                 << " of "
                                 << overlayType
                                 << " Overlay"
                                 << endl;
    TH2F*   hInputHisto             = NULL;
    TH1F*   hOutputMaxHisto         = NULL;
    TH1F*   hOutputMeanHisto        = NULL;
    TH1F*   hOutputMedianHisto      = NULL;
    TH1D*    hTempHisto             = NULL;
    double  max_prop                = 0;
    double  median                  = 0;
    double  mean                    = 0;
    double  max_prop_err            = 0;
    double  median_err              = 0;
    double  mean_err                = 0;

    if (verbosityLevel > 3)
    {
        cout << "...setting pointers to histogramm arrays ";
        cout << " for " << overlayType << "Overlay" << endl;
    }
    if (overlayType.Contains("Edge"))
    {
        hInputHisto        = (CurrentPixel->hEdgeOverlay[pulse_order]);

        //Maximum Propability of Slices
        hOutputMaxHisto    = (CurrentPixel->hPixelEdgeMax[pulse_order]);

        //Mean of Slices
        hOutputMeanHisto   = (CurrentPixel->hPixelEdgeMean[pulse_order]);

        //Median of Slices
        hOutputMedianHisto = (CurrentPixel->hPixelEdgeMedian[pulse_order]);
    }
    else if (overlayType.Contains("Maximum"))
    {
        hInputHisto        = (CurrentPixel->hMaxOverlay[pulse_order]);

        //Maximum Propability of Slices
        hOutputMaxHisto    = (CurrentPixel->hPixelMax[pulse_order]);

        //Mean of Slices
        hOutputMeanHisto   = (CurrentPixel->hPixelMean[pulse_order]);

        //Median of Slices
        hOutputMedianHisto = (CurrentPixel->hPixelMedian[pulse_order]);
    }
    else
    {
        cout << endl << overlayType << "Unknown Overlay Method-->aborting" << endl;
        return;
    }
    if (verbosityLevel > 3)
    {
        cout << "...done " << endl;
    }

    if (verbosityLevel > 2)
    {
        cout << "\t...calculation of "
             << "pulse order " << pulse_order << endl;
    }

//    if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;

    int slice_min  = hInputHisto->GetXaxis()->GetFirst();
    int slice_max  = hInputHisto->GetXaxis()->GetLast();

    if (verbosityLevel > 2)
    {
        cout << "\t...looping over slice range " << slice_min
             << " to " << slice_max << endl;
    }

    for (Int_t slice=slice_min;slice<=slice_max;slice++)
    {

        hTempHisto  = hInputHisto->ProjectionY("",slice,slice,"o");

        //containers for errors
        double slMean[3]    ={0,0,0};
        double slError[3]   ={0,0,0};

        //calculate Errors with bootstrapping
        CalculateErrorsWithBootstrapping(hTempHisto, numBootstrapIt, slMean, slError);

        if (verbosityLevel > 3)
        {
            cout << "\t\t...calculating maxProb of slice " << slice << endl;
        }

        //Get maximum of slice's distribution
//        max_prop    = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
                max_prop    = slMean[0];

        //improve result by< fitting gaussian to slices distribution
//        TF1 gaus("fgaus", "gaus", max_prop-30, max_prop+30);
//        hTempHisto->Fit("fgaus", "QRN0");
//        max_prop        = gaus.GetParameter(1);

        //calculate error of max prop
//        max_prop_err    = gaus.GetParameter(2); //RMS of gaus fit of slices distribution
        max_prop_err    = slError[0];         //error calculated with bootstrapping
//        max_prop_err    = hTempHisto->GetRMS();   //RMS of slice

        if (verbosityLevel > 3)
        {
            cout << "\t\t...calculating Median of slice " << slice << endl;
        }
//        median      = MedianOfH1(hTempHisto);
        median      = slMean[1];
        median_err  = slError[1]; //error from bootstraping

        if (verbosityLevel > 3) cout << "\t\t...calculating Mean of slice " << slice << endl;


//        mean        = hTempHisto->GetMean();
        mean        = slMean[2];
//        mean_err    = hTempHisto->GetRMS()/hTempHisto->GetEntries(); //RMS of slice
        mean_err    = slError[2];         //error from bootstraping

        if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
        hOutputMaxHisto->SetBinContent(slice, max_prop );
        hOutputMaxHisto->SetBinError(    slice,  max_prop_err);

        if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
        hOutputMeanHisto->SetBinContent(slice, mean );
        hOutputMeanHisto->SetBinError(    slice,  mean_err);

        if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
        hOutputMedianHisto->SetBinContent(  slice,  median );
        hOutputMedianHisto->SetBinError(    slice,  median_err);
        delete hTempHisto;

    }//Loop over slices

//    ErrorMedianOfTH2Slices(
//            hInputHisto,
//            hOutputMedianHisto,
//            100,  //numIterations
//            verbosityLevel
//            );
//    hOutputMaxHisto->GetXaxis()->SetLimits(
//                0,
//                300
//                );

//    hOutputMeanHisto->GetXaxis()->SetLimits(
//                0,
//                300
//                );

//    hOutputMedianHisto->GetXaxis()->SetLimits(
//                0,
//                300
//                );

    ShiftStartOfHistoInXTo(
            hOutputMaxHisto,
            0
            );

    ShiftStartOfHistoInXTo(
            hOutputMeanHisto,
            0
            );

    ShiftStartOfHistoInXTo(
            hOutputMedianHisto,
            0
            );

    ShiftStartOfHistoInXTo(
            hInputHisto,
            0
            );

}
// end of PlotMaxPropabilityPulse
//----------------------------------------------------------------------------

bool
WritePixelTemplateToCsv(
        Pixel*          CurrentPixel,
        TString         path,
        TString         overlayMethod,
        int             order,
        int             verbosityLevel
        )
{
//    TSystem this_system();
//    this_system.cd(path);
//    cout << this_system.pwd();
//    this_system.mkdir("CSV");
    path = CurrentPixel->CsvFileName( path, overlayMethod, order);

    TH1F*  Max_histo    = NULL;
    TH1F*  Median_histo    = NULL;
    TH1F*  Mean_histo    = NULL;

    if (overlayMethod.Contains("Maximum"))
    {
        Max_histo = CurrentPixel->hPixelMax[order];
            Median_histo = CurrentPixel->hPixelMedian[order];
            Mean_histo = CurrentPixel->hPixelMean[order];
    }
    else if (overlayMethod.Contains("Edge"))
    {
        Max_histo = CurrentPixel->hPixelMax[order];
        Median_histo = CurrentPixel->hPixelMedian[order];
        Mean_histo = CurrentPixel->hPixelMean[order];
    }
    else
    {
        cout << endl << "Unknown Overlay Method-->aborting" << endl;
        return 1;
    }


    Int_t nbins = Max_histo->GetXaxis()->GetNbins();

    if (verbosityLevel > 2)
    {
        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: "
        << overlayMethod << endl;
    out << "### Slice's Amplitude determined by calculating the " << endl
        << "### value of maximum propability of slice -> AmplitudeMax " << endl
        << "### mean of slice -> AmplitudeMean " << endl
        << "### median of slice -> AmplitudeMedian " << endl
        << "### for each slice" << endl;
    out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl
        << endl;

    out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;

    for (int slice=1;slice<=nbins;slice++)
    {
        out << slice << "," ;
        out << Max_histo->GetBinContent(slice) << ",";
        out << Mean_histo->GetBinContent(slice) << ",";
        out << Median_histo->GetBinContent(slice) << endl;
    }
    out.close();
    if (verbosityLevel > 2) cout << "...csv file closed" << endl;
    return 0;
}

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
//FitFallingEdge(
//        TString         name,
//        TH1F*           histo,
//        double          xMin,
//        double          xMax,
//        double*         parameters
//        )
//{
//    TF1* polExpFit = new TF1(name, PolExp, xMin, xMax, 3 );
//    polExpFit->SetParNames("Pol0", "Slope", "Shift");
//    polExpFit->SetLineColor(kGreen);
//    histo->Fit(polExpFit, "+RWM");
//    polExpFit->GetParameters(parameters);
//}

//void
//FitRisingEdge(
//        TString         name,
//        TH1F*           histo,
//        double          xMin,
//        double          xMax,
//        double*         parameters
//        )
//{
//    TF1* polExpFit = new TF1(name, NegPolExp, xMin, xMax, 3 );
//    polExpFit->SetParNames("Pol0", "Slope", "Shift");
//    polExpFit->SetLineColor(kRed);
//    histo->Fit(polExpFit, "+RWWM");
//    polExpFit->GetParameters(parameters);
//}

//double
//NegPolExp(
//        double*         x,
//        double*         par
//        )
//{
//    return par[0]+(-1)*TMath::Exp(par[1]+par[2]*x[0]);
//}

//double
//PolExp(
//        double*         x,
//        double*         par
//        )
//{
//    return
////            par[0]+
//            TMath::Exp(par[1]+par[2]*x[0]);
//}

//double
//ChargeDiode(
//        double         time,
//        double         chargeVoltage,
//        double         impedance,
//        double         capacity
//        )
//{
//    return chargeVoltage*(1 - TMath::Exp(time/(impedance*capacity)));
//}

//double
//UnChargeDiode(
//        double*         time,
//        double*         chargeVoltage,
//        double*         timeConstant
//        )
//{
//   return chargeVoltage[0]+TMath::Exp(chargeVoltage[1]+timeConstant[2]*time[0]);
////    return chargeVoltage[0] * (TMath::Exp(time[0]/timeConstant[0]));
//}

//double
//template_function(
//        double*         input_x,
//        double*         par)
//{
//    double returnval = 0.0;

//    // I introduce a few names
//   // double shift = par[0];
//    double bsl = par[0];
//    double beginOfRisingEdge = par[1];
//    double p0RisingEdge = par[6];
//    double p1RisingEdge = par[7];
//    double p2RisingEdge = par[8];
//    double p3RisingEdge = par[9];
//    double endOfRisingEdge = par[2];
////    double pOFallingEdge = par[3];
////    double expPar1FallingEdge = par[4];
////    double expPar1FallingEdge = par[5];
//    /*
//    bool couted_once = false;
//    if not couted_once
//    {
//        couted_once = true;
//        cout << "shift:" << shift << endl;
//        cout << "bsl:" << bsl << endl;
//        cout << "expars:" << endl;
//        cout << "\t factor:" << exppar[0] << endl;
//        cout << "\t tau:" << exppar[1] << endl;
//        cout << "\t t0:" << exppar[2] << endl;
//        cout << "pol3pars:" << endl;
//        cout << "p[0] + x p[1] + x^2 p[2] + x^3 p[3]" << endl;
//        cout << pol3par[0] << "\t" << pol3par[1] << "\t" << pol3par[2] << "\t" << pol3par[3] << endl;
//        cout << "ranges:" << endl;
//        cout << "begin of pol3: " << range[0] << endl;
//        cout << "begin of exp:  " << range[1] << endl;
//    }
//    */
//    double x = input_x[0];

//    // the baseline is added everywhere.
//    returnval += bsl;

//    if ( (x > beginOfRisingEdge) && (x <= endOfRisingEdge) )
//    {
//        // from this point on the pol3 is added
//        returnval += p0RisingEdge;
//        returnval += p1RisingEdge * x;
//        returnval += p2RisingEdge * x*x;
//        returnval += p3RisingEdge * x*x*x;
//    }
//    else if ( x > endOfRisingEdge )
//    {
//        // from this point on the exp-func is added
////        returnval += exppar[0] * TMath::Exp( exppar[1] * ( x - exppar[2] ) );
//        returnval += PolExp(input_x, par+3);
//    }

//    return returnval;
//}


void
ShiftHistoInX(
        TH1*                histo,
        float               shift
        )
{
    histo->GetXaxis()->SetLimits(
        histo->GetXaxis()->GetBinLowEdge( histo->GetXaxis()->GetFirst() ) + shift,
        histo->GetXaxis()->GetBinUpEdge(  histo->GetXaxis()->GetLast() )  + shift
                              );
}

void
ShiftStartOfHistoInXTo(
        TH1*                histo,
        float               value
        )
{
    float min = histo->GetXaxis()->GetBinLowEdge( histo->GetXaxis()->GetFirst() );
    float shift = value - min;

    ShiftHistoInX( histo, shift);
}