#include <iostream>
#include <fstream>

#include "templateextractors.h"
#include <stdlib.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   max_value_of_slice  = 0;
    TH1D*   hTempHisto          = NULL;
    int     nbins               = inputHisto->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
        hTempHisto          = inputHisto->ProjectionY( "", TimeSlice,TimeSlice);
        max_value_of_slice  = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );

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

        outputHisto->SetBinContent(TimeSlice, max_value_of_slice );
    }

    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(kRed);
    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, "RWM");
    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*         impedance,
        double*         capacity,
        )
{
    return chargeVoltage*(TMath::Exp(time/(impedance*capacity)));
}

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

double
PulseFunction(
        double*         time,
        double*         baseline,
        double*         risingChargeVoltage,
        double*         risingImpedance,
        double*         risingCapacity,
        double*         fallingChargeVoltage,
        double*         fallingImpedance,
        double*         fallingCapacity
        )
{
    double returnValue  = 0.0;
    returnValue         += baseline;
    returnValue         +=  ChargeDiode(
                                    time,
                                    risingChargeVoltage,
                                    risingImpedance,
                                    risingCapacity);
    returnValue         +=  UnChargeDiode(
                                    time,
                                    fallingChargeVoltage,
                                    fallingImpedance,
                                    fallingCapacity);
    return returnValue;
}

void
FitPulse(
        TString         name,
        TH1F*           histo,
        double          xMin,
        double          xMax,
        double*         parameters
        )
{
    TF1* pulseFunction = new TF1(name, PulseFunction, xMin, xMax, 7 );
    pulseFunction->SetParNames(
                "Baseline",
                "Charge-Voltage of rising Edge",
                "Impedance for rising Edge",
                "Capacity for rising Edge",
                "Charge-Voltage of falling Edge",
                "Impedance for falling Edge",
                "Capacity for falling Edge"
                );
    pulseFunction->SetLineColor(kRed);
    histo->Fit(pulseFunction, "RWM");
    pulseFunction->GetParameters(parameters);
    return 0;
}