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