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

#include "templateextractors.h"

using namespace std;


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 timeslices 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_t MedianOfH1 (
        TH1*            inputHisto
        )
{
   //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(i+1);
      y[i] = inputHisto->GetBinContent(i+1);
   }
   Double_t median = TMath::Median(nbins,x,y);
   delete [] x;
   delete [] y;
   return median;
}
// end of PlotMedianEachSliceOfPulse
//----------------------------------------------------------------------------

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;

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

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

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

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

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

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

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 TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {

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

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

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

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

void
ExtractPulseTemplate(
        Pixel*          CurrentPixel,
        TString         overlayType,
        int             pulse_order,
        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;
    TH1*    hTempHisto              = NULL;
    float   max_prop                = 0;
    float   median                  = 0;
    float   mean                    = 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;

    for (Int_t TimeSlice=1;TimeSlice<=300;TimeSlice++)
    {

        hTempHisto  = hInputHisto->ProjectionY("",TimeSlice,TimeSlice);

        if (verbosityLevel > 3)
        {
            cout << "\t\t...calculating maxProb of slice " << TimeSlice << endl;
        }
        max_prop    = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );

        if (verbosityLevel > 3)
        {
            cout << "\t\t...calculating Median of slice " << TimeSlice << endl;
        }
        median      = MedianOfH1(hTempHisto);

        if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << TimeSlice << endl;
        mean        = hTempHisto->GetMean();

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

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

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

    }//Loop over Timeslices
}
// 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 TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    {
        out << TimeSlice << "," ;
        out << Max_histo->GetBinContent(TimeSlice) << ",";
        out << Mean_histo->GetBinContent(TimeSlice) << ",";
        out << Median_histo->GetBinContent(TimeSlice) << 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;
//}