#include "TCanvas.h"
#include "TH1.h"
#include "TH1F.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
#include "TVirtualFFT.h"
#include <fstream>
#include <iostream>
#include "TStyle.h"





void UpSampleTrace(){
  
  cout<<"Up-sampling trace"<<endl;
   //read in trace
  gStyle->SetOptStat(0);
  double *x = new double[1000];
  float *y = new float[1000];
  int iter = 0, TankId;
  long int EventId;
  double Signal, Distance, RiseTime, Energy, EnergyError, Theta, Zeta;
  double S1000, S1000Error;
  double ConvolvedRiseTime, DeConvolvedRiseTime;
  double ConvolvedSignal, DeconvolvedSignal;
  
  ifstream Data;
  Data.open("fft_forum.txt");
 
  while (Data >> x[iter] >> y[iter] >> Signal >> EventId >> TankId >> Distance 
	 >> RiseTime >> Zeta >> Energy >> EnergyError >> Theta 
	 >>S1000 >> S1000Error)
    iter++;
  
  //load it into a histogram and draw
  TCanvas* NewCanvas=new TCanvas(Form("%d %d", EventId, TankId),Form("%d %d", EventId, TankId)); 
  NewCanvas->Divide(2,1);
  NewCanvas->cd(1);
  TH1F *h = new TH1F(Form("h %d %d", EventId, TankId),"FADC trace",iter,x[0],x[iter]);
  //h->GetXaxis()->SetRange(x[0], x[iter]);
  for (int i=0;i<iter;i++)
    h->SetBinContent(i+1,y[i]/3);
  h->GetXaxis()->SetTitle("t [ns]");
  h->GetYaxis()->SetTitle("Signal [VEM]");
  h->Draw();
  
  
  //transform into frquency domain
  NewCanvas->cd(2);
  TH1F *hTrans = new TH1F("h in frequecy space", "h in frequecy space", iter, 0, iter);
  hTrans = (TH1F*)h->FFT(hTrans, "MAG R2C M");
  hTrans->Draw();
  hTrans->SetStats(kFALSE);
  

  Double_t *output = new Double_t[2*((iter)/2+1)];

  TVirtualFFT *current = TVirtualFFT::GetCurrentTransform();
  current->GetPoints(output);

  TVirtualFFT *inverse = TVirtualFFT::FFT(1, &iter, "C2R M K");
  inverse->SetPoints(output);
  inverse->Transform();
  inverse->GetPoints(output);
  for (Int_t i=0; i<10; i++){
     printf("output[%d]=%f, h[%d]=%f\n", i, output[i]/Double_t(iter), i+1, h->GetBinContent(i+1));
  }


  delete [] output;
  delete current;
  delete inverse;
  //for (Int_t i=0; i< 
  //transform back into time domain - should recovere same trace
//   TH1F *hInvTrans = new TH1F("h in time", "h in time", iter, x[0],x[iter] );
//   hInvTrans = hTrans->FFT(hInvTrans, "MAG C2R M");
//   hInvTrans->Draw();
//   hInvTrans->SetStats(kFALSE);
}