#include <math.h>
#include <iostream>
#include "Filter.h"

ClassImp(Filter);
ClassImp(Filter::PeakIterator);

#ifndef SQR_FUNC
#define SQR_FUNC 1
Double_t sqr(Double_t x) { return x*x; }
#endif

Filter::Filter(const TH1D &h1)
 : fSize(0), fSeries(0)
{
   fLength = h1.GetXaxis()->GetNbins();
   fCoef = new Double_t[fLength];
   for (Int_t i=0; i < fLength; ++i) {
      fCoef[i] = h1.GetBinContent(i+1);
   }
}

Filter::Filter(const Filter &filt)
 : fSize(0), fSeries(0)
{
   fLength = filt.fLength;
   for (Int_t i=0; i < fLength; ++i) {
      fCoef[i] = filt.fCoef[i];
   }
}

Filter::Filter(Int_t length)
 : fSize(0), fSeries(0)
{
   Int_t nPed = 4;
   Int_t nFlat = 2;
   Double_t t0 = 1.5;
   Double_t t1 = 5;
   fLength = length;
   fCoef = new Double_t[length];

   Double_t sum_pre=0;
   Double_t sum_post=0;
   Int_t i=0;
   for (; i<nPed; i++) {
     fCoef[i] = exp(i/t0);
      sum_pre += fCoef[i];
   }
   for (i=nPed; i<nPed+nFlat; i++) {
      fCoef[i] = 1;
      sum_post += fCoef[i];
   }
   for (Int_t t=0; i<length; i++,t++) {
      fCoef[i] = exp(-t/t1);
      sum_post += fCoef[i];
   }
   for (i=0; i<nPed; i++) {
      fCoef[i] *= -sum_post/sum_pre;
   }
}

Filter::~Filter()
{
   if (fSize && fSeries) {
      delete [] fSeries;
   }
}

TH1D* Filter::GetFilterHist()
{
   TH1D *h = new TH1D("filter","linear filter",fLength,0,fLength);
   for (Int_t i=0; i < fLength; ++i) {
      h->SetBinContent(i+1,fCoef[i]);
   }
   return h;
}

Int_t Filter::GetFilterData(Double_t *array, Int_t length)
{
   Int_t limit = (length < fLength)? length : fLength;
   for (Int_t i=0; i<limit; ++i) {
      array[i] = fCoef[i];
   }
   return limit;
}

void Filter::SetEdgeFilter(Int_t len)
{
   len = (len == 0)? fLength : len;
   if (len != fLength) {
      delete [] fCoef;
      fCoef = new Double_t[len];
      fLength = len;
   }
   Double_t sum=0;
   Int_t mid=len/2;
   Double_t width=mid/2;
   for (Int_t i=0; i < mid; ++i) {
      fCoef[i] = fCoef[len-i-1] = exp((i-mid+1)/width);
      sum += fCoef[i]*((len-i-1 > mid)? 2 : 1);
   }
   fCoef[mid] = -sum;
}

TH1D *Filter::GetTransformHist()
{
   if (fSize == 0) {
      return 0;
   }
   TH1D *h = new TH1D("transform","filtered data series",fSize,0,fSize);
   for (Int_t i=0; i < fSize; ++i) {
      h->SetBinContent(i+1,fSeries[i]);
   }
   return h;
}

Int_t Filter::GetTransformData(Double_t *array, Int_t length)
{
   Int_t limit = (length < fSize)? length : fSize;
   if (limit == 0) {
      return 0;
   }
   for (Int_t i=0; i < limit; ++i) {
      array[i] = fSeries[i];
   }
   return limit;
}

void Filter::Transform(const TH1 *hin)
{
   Int_t length = hin->GetXaxis()->GetNbins();
   Double_t *x = new Double_t[length];
   for (Int_t i=0; i < length; ++i) {
      x[i] = hin->GetBinContent(i+1);
   }
   Transform(x,length);
   delete [] x;
}

void Filter::Transform(const Double_t *array, Int_t length)
{
   if (fSize && fSeries) {
      delete [] fSeries;
   }
   fSize = length-fLength+1;
   fSize = (fSize < 0)? 0 : fSize;
   if (fSize == 0) {
      return;
   }
   fSeries = new Double_t[fSize];
   for (Int_t i0=0; i0 < fSize; ++i0) {
      fSeries[i0] = 0;
      for (Int_t i=0; i < fLength; ++i) {
         fSeries[i0] += array[i0+i]*fCoef[i];
      }
   }
}

void Filter::Transform(const Float_t *array, Int_t length)
{
   if (fSize && fSeries) {
      delete [] fSeries;
   }
   fSize = length-fLength+1;
   fSize = (fSize < 0)? 0 : fSize;
   if (fSize == 0) {
      return;
   }
   fSeries = new Double_t[fSize];
   for (Int_t i0=0; i0 < fSize; ++i0) {
      fSeries[i0] = 0;
      for (Int_t i=0; i < fLength; ++i) {
         fSeries[i0] += array[i0+i]*fCoef[i];
      }
   }
}

void Filter::Transform(const Int_t *array, Int_t length)
{
   if (fSize && fSeries) {
      delete [] fSeries;
   }
   fSize = length-fLength+1;
   fSize = (fSize < 0)? 0 : fSize;
   if (fSize == 0) {
      return;
   }
   fSeries = new Double_t[fSize];
   for (Int_t i0=0; i0 < fSize; ++i0) {
      fSeries[i0] = 0;
      for (Int_t i=0; i < fLength; ++i) {
         fSeries[i0] += array[i0+i]*fCoef[i];
      }
   }
}

void Filter::Transform(const Short_t *array, Int_t length)
{
   if (fSize && fSeries) {
      delete [] fSeries;
   }
   fSize = length-fLength+1;
   fSize = (fSize < 0)? 0 : fSize;
   if (fSize == 0) {
      return;
   }
   fSeries = new Double_t[fSize];
   for (Int_t i0=0; i0 < fSize; ++i0) {
      fSeries[i0] = 0;
      for (Int_t i=0; i < fLength; ++i) {
         fSeries[i0] += array[i0+i]*fCoef[i];
      }
   }
}

void Filter::Streamer(TBuffer &b)
{
   if (b.IsReading()) {
      TNamed::Streamer(b);
      b >> fLength;
      for (Int_t i=0; i < fLength; ++i) {
         b >> fCoef[i];
      }
      b >> fSize;
      for (Int_t i=0; i < fSize; ++i) {
         b >> fSeries[i];
      }
   }
   else if (b.IsWriting()) {
      TNamed::Streamer(b);
      b << fLength;
      for (Int_t i=0; i < fLength; ++i) {
         b << fCoef[i];
      }
      b << fSize;
      for (Int_t i=0; i < fSize; ++i) {
         b << fSeries[i];
      }
   }
}

Filter::PeakIterator &Filter::GetPeakIterator(Double_t thresh)
{
   PeakIterator *piter = new PeakIterator(this,thresh);
   return *piter;
}

Filter::PeakIterator::PeakIterator(Filter *filt, Double_t thresh)
{
   fThresh = thresh;
   fFilter = filt;
   fPeak.pos = -1;
}

Filter::PeakIterator::PeakIterator(const PeakIterator &iter)
{
   fThresh = iter.fThresh;
   fFilter = iter.fFilter;
   fPeak.pos = iter.fPeak.pos;
   fPeak.mag = iter.fPeak.mag;
}

Filter::PeakIterator::~PeakIterator()
{
}

void Filter::PeakIterator::Reset()
{
   fPeak.pos = -1;
}

Filter::PeakIterator Filter::PeakIterator::Begin()
{
   PeakIterator biter(*this);
   biter.fPeak.pos = -1;
   return biter;
}

Filter::PeakIterator Filter::PeakIterator::End()
{
   PeakIterator diter(*this);
   diter.fPeak.pos = fFilter->fSize;
   return diter;
}

Filter::PeakIterator &Filter::PeakIterator::operator=(PeakIterator &iter)
{
   fThresh = iter.fThresh;
   fFilter = iter.fFilter;
   fPeak.pos = iter.fPeak.pos;
   fPeak.mag = iter.fPeak.mag;
   return *this;
}

Filter::Peak_t &Filter::PeakIterator::operator*()
{
   return fPeak;
}

Bool_t Filter::PeakIterator::operator==(Filter::PeakIterator &iter)
{
   if ((fThresh != iter.fThresh) ||
       (fFilter != fFilter) ||
       (fPeak.pos != iter.fPeak.pos)) {
      return false;
   }
   else {
      return true;
   }
}

Bool_t Filter::PeakIterator::operator!=(Filter::PeakIterator &iter)
{
   if (*this == iter) {
      return false;
   }
   else {
      return true;
   }
}

Filter::PeakIterator &Filter::PeakIterator::operator++()
{
   Int_t ipeak = -1;
   Int_t i = fPeak.pos;
   Double_t peak = fThresh;
   while (++i < fFilter->fSize) {
      Double_t samp = fFilter->fSeries[i];
      if (samp >= peak) {
         peak = samp;
         ipeak = i;
      }
      else if (samp < 0 && ipeak != -1) {
         break;
      }
   }
   if (ipeak == -1) {
      fPeak.mag = 0;
      fPeak.pos = -1;
      fPeak.zero = -1;
   }
   else {
      fPeak.mag = peak;
      fPeak.pos = i;
      Double_t samp_neg = fFilter->fSeries[i];
      Double_t samp_pos = fFilter->fSeries[i-1];
      fPeak.zero = i+samp_neg/(samp_pos-samp_neg);
   }
   return *this;
}

Filter::PeakIterator Filter::PeakIterator::operator++(int)
{
   PeakIterator me(*this);
   ++(*this);
   return me;
}

void Filter::PeakIterator::Streamer(TBuffer &b)
{
   std::cerr << "Error - cannot call Streamer method for PeakIterator!"
             << std::endl;
}