#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,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of max. propability"
                                 << endl;
    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(
                    CurrentPixel->hMaxOverlay[pulse_order],
                    CurrentPixel->hPixelMax[pulse_order],
                    verbosityLevel);

        CalcMaxPropabilityOfSlice(
                    CurrentPixel->hEdgeOverlay[pulse_order],
                    CurrentPixel->hPixelEdgeMax[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
CalcMedianOfSlice(
        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("Max"))
    {
        hInputHistoArray    = CurrentPixel->hMaxOverlay;
        hOutputHistoArray   = CurrentPixel->hPixelMedian;
    }

    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
//----------------------------------------------------------------------------

bool
WritePixelTemplateToCsv(
        Pixel*          CurrentPixel,
        TString         path,
        TString         overlayMethod,
        int             order,
        int             verbosityLevel
        )
{
    path = CurrentPixel->CsvFileName( path, overlayMethod, order);

    TH1F**  Max_histo_array    = NULL;
    TH1F**  Median_histo_array    = NULL;
    TH1F**  Mean_histo_array    = NULL;

    if (overlayMethod.Contains("Maximum"))
    {
            Max_histo_array = CurrentPixel->hPixelMax;
            Median_histo_array = CurrentPixel->hPixelMedian;
            Mean_histo_array = CurrentPixel->hPixelMean;
    }
    else if (overlayMethod.Contains("Edge"))
    {
        Max_histo_array = CurrentPixel->hPixelMax;
        Median_histo_array = CurrentPixel->hPixelMedian;
        Mean_histo_array = CurrentPixel->hPixelMean;
    }
    else
    {
        cout << endl << "Unknown Overlay Method-->aborting" << endl;
        return 1;
    }


    Int_t nbins = Max_histo_array[order]->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: "
        << 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_array[order]->GetBinContent(TimeSlice) << ",";
        out << Mean_histo_array[order]->GetBinContent(TimeSlice) << ",";
        out << Median_histo_array[order]->GetBinContent(TimeSlice) << endl;
    }
    out.close();
    if (verbosityLevel > 0) cout << "...file closed" << endl;
    return 0;
}