//////////////////////////////////////////////////////////////
// Functions for FPulsetemplate to calculate
// values for slices of TH2 histograms
//////////////////////////////////////////////////////////////

#include "slicecalculation.h"
#include <TROOT.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>

//compute the median for 1-d histogram h1
Double_t
MedianOfH1(
        TH1*            h1,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 4) cout << endl
                                 << "...calculating median of given TH1"
                                 << endl;
   Int_t nbins = h1->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] = h1->GetXaxis()->GetBinCenter(i+1);
      y[i] = h1->GetBinContent(i+1);
   }
   Double_t median = TMath::Median(nbins,x,y);

   delete [] x;
   delete [] y;

   if (verbosityLevel > 4) cout << "...done" << endl;
   return median;
}

//compute the median for each slice of a given TH2 histogram
void
PlotMedianEachSliceOfPulse(
        TH2*            phInputHistogram,
        TH1*            phOutputHistogram,
        int             fitdata,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of slice's Median"
                                 << endl;

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

    for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    {
       TH1 *hProjection = phInputHistogram->ProjectionY("",TimeSlice,TimeSlice);
       double median = MedianOfH1(hProjection);

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

       delete hProjection;

       phOutputHistogram->SetBinContent(TimeSlice, median );
//       hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
    }

    if (verbosityLevel > 2) cout << "\t...done" << endl;

    if (fitdata)
    {
       FitMaxPropabilityPuls(
                phOutputHistogram,
                verbosityLevel
                );
    }
}

//compute the mean for each slice of a given TH2 histogram
void
PlotMeanEachSliceOfPulse(
        TH2*            phInputHistogram,
        TH1*            phOutputHistogram,
        int             fitdata,
        int             verbosityLevel
        )
{
    if (verbosityLevel > 2) cout << endl
                                 << "...calculating pulse shape of slice's mean"
                                 << endl;

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

    for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    {
       TH1 *hProjection = phInputHistogram->ProjectionY("",TimeSlice,TimeSlice);
       double mean = hProjection->GetMean();

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

       delete hProjection;

       phOutputHistogram->SetBinContent(TimeSlice, mean );
//       hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
    }

    if (verbosityLevel > 2) cout << "\t...done" << endl;

    if (fitdata)
    {
       FitMaxPropabilityPuls(
                phOutputHistogram,
                verbosityLevel
                );
    }
}

//fit a function to the given histogram
void
FitMaxPropabilityPuls(
        TH1*           phInputHistogram,
        float           left_border,
        float           right_border,
        int             verbosityLevel
        )
    {
      if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
      if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
      phInputHistogram->Fit("landau", "", "", left_border, right_border);
      if (verbosityLevel > 2) cout << "...done" << endl;
    }