// Author: Andrea Palaia  24/08/09

#include <stdlib.h>

#include "Riostream.h"
#include "TVectorD.h"
#include "TGraph.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TMath.h"
#include "Waveform.h"

Int_t index(TVectorD v, Double_t val) {

  Int_t index = 0;
  Int_t size = v.GetNoElements();

  for (Int_t i=0; i<size; i++){
    if(v(i)==val) {
      index = i;
      break;
    }
  }
  return index;
}

Int_t indexClosestVal(TVectorD v, Double_t val, int start = 0, int stop = 0) {

  Double_t diff = fabs(v(0)-val);
  Int_t index = 0;

  if(stop==0) stop = v.GetNoElements();

  for (int i=start; i<stop; i++){
    if (fabs(v(i)-val)<diff) {
      diff=fabs(fabs(i)-val);
      index = i;
    }
  }
  return index;
}

ClassImp(Waveform)

//______________________________________________________________________________
  Waveform::Waveform(){
  // Waveform default constructor.
}

//______________________________________________________________________________
Waveform::Waveform(Double_t v[], UInt_t size) {
  SetWaveform(v, size);
}

//______________________________________________________________________________
Waveform::Waveform(Float_t v[], UInt_t size) {
  SetWaveform(v, size);
}

//______________________________________________________________________________
Waveform::Waveform(TVectorD v) {
  SetWaveform(v);
}

//______________________________________________________________________________
void Waveform::SetWaveform(Double_t v[], UInt_t size) {
  wave.ResizeTo(size);
  for(UInt_t i=0;i<size;i++) {
    wave(i) = v[i];
  }
  SetSampleFreq(1);
}

//______________________________________________________________________________
void Waveform::SetWaveform(Float_t v[], UInt_t size) {
  wave.ResizeTo(size);
  for(UInt_t i=0;i<size;i++) {
    wave(i) = (Double_t)v[i];
  }
  SetSampleFreq(1);
}

//______________________________________________________________________________
void Waveform::SetWaveform(TVectorD v) {
  wave.ResizeTo(v.GetNoElements());
  wave.Add(v);
  SetSampleFreq(1);
}

//______________________________________________________________________________
TVectorD Waveform::GetWaveform() {
  return wave;
}

//______________________________________________________________________________
void Waveform::SetSampleFreq(Double_t f) {
  // Set the samplig frequency during acquisition: this is important in order to have the proper time scale (plots and fwhm).

  sampleFreq = f;
}

//______________________________________________________________________________
Double_t Waveform::GetSampleFreq() {
  return sampleFreq;
}

//______________________________________________________________________________
void Waveform::SetLabel(TString str) {
  label = str;
}

//______________________________________________________________________________
TString Waveform::GetLabel() {
  return label;
}

//______________________________________________________________________________
void Waveform::SetOffset() {
  // This calculates the zero level of the waveform and subract it to the spectrum (the old spectrum will be overwritten).
  // The zero level is calculated averaging a number of points at the beginning of the spectrum (the first point is skipped).
  // If the spectrum is recognised to be noise the first 50 points will be averaged, otherwise the first 10
  // (this is to prevent averaging over part of the signal).

  UInt_t nPoints = 10;
  if(IsNoise()) nPoints = 50;
  Double_t sum = 0;
  for(UInt_t j=1;j<nPoints;j++) {
    sum += wave(j);
  }
  offset = sum/nPoints;
  for(UInt_t j=0;j<abs(wave.GetNoElements());j++) {
    wave(j) = wave(j)-offset;
  }
}

//______________________________________________________________________________
Double_t Waveform::GetOffset() {
  return offset;
}

//______________________________________________________________________________
void Waveform::SetUnixtime(ULong_t t) {
  unixtime = t;
}

//______________________________________________________________________________
ULong_t Waveform::GetUnixtime() {
  return unixtime;
}

//______________________________________________________________________________
TVectorD Waveform::Smooth() {
  // This return a TVectorD which contains a smoothed spectrum.
  Int_t size = wave.GetNoElements();
  TVectorD smoothed_wave(size);
  for(Int_t j=0;j<size;j++) {
    smoothed_wave(j) = 0.0;
  }
  for(Int_t j=0;j<size-2;j++) {
    smoothed_wave(j+1) = (wave(j)+wave(j+1)+wave(j+2))/3;
  }
  return smoothed_wave;
}

//______________________________________________________________________________
Double_t Waveform::Fwhm() {
  // Calculates the Full Width half Maximun of the spectrum.
  return (GetIndexAt(0.5,1) - GetIndexAt(0.5))*(1/sampleFreq);

}

//______________________________________________________________________________
void Waveform::Plot(Option_t* option, Color_t color) {
  TString opt = option;
  Int_t size = wave.GetNoElements();
  TVectorD time(size);
  for(Int_t i=0;i<size;i++) {
    time(i)=i/sampleFreq;
  }
  if(opt=="") opt.Append("AL");
  if(opt.Contains("A")) TCanvas* c = new TCanvas();
  TGraph* gr1 = new TGraph(time, wave);
  gr1->SetLineColor(color);
  gr1->SetLineWidth(2);
  gr1->GetXaxis()->SetTitle("s");
  gr1->Draw(opt.Data());
}

//______________________________________________________________________________
Double_t Waveform::GetMax() {
  // Gets the maximum value of the spectrum.

  return wave.Max();
}

//______________________________________________________________________________
Double_t Waveform::Integral() {
  Double_t integral = 0;
  Int_t size = wave.GetNoElements();
  for(Int_t j=0;j<size-1;j++) {
    integral += wave(j) + (fabs(wave(j)-wave(j+1))/2);
  }
  return integral*(1/sampleFreq);
}

//______________________________________________________________________________
Int_t Waveform::GetIndexAt(Double_t frac, Bool_t falling_edge) {
  Int_t index = -1;
  Double_t val = GetMax() * frac;
  Int_t start, stop;
  if(falling_edge) {
    start = IndexOf(GetMax());
    stop = wave.GetNoElements()-1;
  } else {
    start = 0;
    stop = IndexOf(GetMax());
  }
  for(Int_t j=start;j<stop;j++) {
    if(!falling_edge&&wave(j)<=val&&wave(j+1)>val) index = j;
    if(falling_edge&&wave(j)>=val&&wave(j+1)<val) index = j;
  }
  return index;
}

//______________________________________________________________________________
Double_t Waveform::GetRisetime() {
  return (GetIndexAt(0.9) - GetIndexAt(0.1))*(1/sampleFreq);
}

//______________________________________________________________________________
Double_t Waveform::GetFalltime() {
  return (GetIndexAt(0.1,1) - GetIndexAt(0.9,1))*(1/sampleFreq);
}

//______________________________________________________________________________
Double_t Waveform::GetMean() {
  Double_t mean = 0;
  Int_t start = GetIndexAt(0.5);
  Int_t stop = GetIndexAt(0.5,1);
  if(start==-1||stop==-1) mean = 0;
  else {
    for(Int_t j=start;j<stop;j++) {
      mean+=wave(j);
    }
  }
  return mean/(stop-start);
}

//______________________________________________________________________________
Double_t Waveform::GetRMS() {
  Double_t ssum = 0;
  Int_t start = GetIndexAt(0.5);
  Int_t stop = GetIndexAt(0.5,1);
  for(Int_t j=start;j<stop;j++) {
    ssum+=pow(wave(j),2);
  }
  return sqrt(ssum/(stop-start));
}

//______________________________________________________________________________
Int_t Waveform::IndexOf(Double_t val) {
  Int_t index = -1;
  for(Int_t j=0;j<wave.GetNoElements();j++) {
    if(wave(j)==val) index = j;
  }
  return index;
  }

//______________________________________________________________________________
Bool_t Waveform::IsSaturated() {
  //Check whether the spectrum is constant for a certain number of points.

  //number of points to be consecutive for the waveform to be said that it is satured
  UInt_t flatPoints = 10;
  Bool_t flag = false;
  UInt_t counter = 0;
  for(UInt_t j=1;j<wave.GetNoElements()-flatPoints;j++) {
    for(UInt_t k=0;k<flatPoints;k++) {
      if(wave(j)==wave(j+k)) counter++;
    }
    if(counter>10) {
      flag = true;
      break;
    } else counter = 0;
  }
  return flag;
}

//______________________________________________________________________________
Bool_t Waveform::IsNoise() {
  Bool_t flag = false;
  if(Integral()*GetSampleFreq()<1) flag = true;
  return flag;
}

  //______________________________________________________________________________
void Waveform::Shift(Int_t points) {

}