int main(){


  float *x = new float[1000];
  float *y = new float[1000];//clenaed trace
  float *z = new float[1000];//un-cleaned trace
  float *yc = new float[1000];//clenaed and deconvolved trace
  int iter = 0, TankId, IsMisfired;
  long int EventId;
  double Signal,  Distance, RiseTime, Energy, EnergyError, Theta, Zeta;
  double S1000, S1000Error, DeconvolvedSignal, CleanedSignal;
  double T1, T2, T5, T10, T15, T20, T25, T30, T35, T40, T45, T50, T60;
  ifstream Data;
  Data.open("1585422_511.txt");
 
  while (Data >> x[iter] >> y[iter]>> z[iter] >> Signal 
	 >>EventId >> TankId >> Distance 
	 >> RiseTime >> Zeta >> Energy >> EnergyError >> Theta 
	 >> S1000 >> S1000Error){ 
    y[iter] /= 3;
    z[iter] /= 3;
    iter++;
  }
  const int NoBins = iter;
   
  if(Signal > 0){

    TH1F *h = new TH1F(Form("h %d %d", EventId, TankId),"FADC trace",NoBins-1,x[0],x[iter-1]);
    h->GetXaxis()->SetRange(x[0], x[iter]);
  
    for (int i=0;i<NoBins;i++) 
      h->SetBinContent(i+1,y[i]);
    h->GetXaxis()->SetTitle("t [ns]"); 
    h->GetYaxis()->SetTitle("Signal [VEM]");
    
    yc = TraceAnalysis(y, NoBins);//This needs to return the array y to yc  which is the deconvolved one 
    
    TH1F *h2 = new TH1F("h2","",NoBins-1, x[0],x[iter-1]);
    h2->SetLineColor(kBlue);
    for (i=0;i<NoBins;i++)
      h2->SetBinContent(i+1,yc[i]);

    h->Draw(); 
    h2->Draw("same");
    
  }
  
  delete [] x;
  delete [] y;
  delete [] z;
  delete [] yc;

  return 0;


};



float* TraceAnalysis(float y[1000], const int NoBins){
  //Use TSpectrum to find the peak candidates
  TSpectrum *spec = new TSpectrum(1000);
  //Int_t nfound = spec->Search(h,2,"");
  
  //set SPR in an array;
  float *SPR = new float[NoBins];//define the response function of the system
  for (int j = 0; j < NoBins; ++j)
    //this is the response function
    SPR[j]=1*(exp(-1*((j*25)+12.5)/67.0) - exp(-1*((j*25)+12.5)/13.0));
  
  //deconvolute SPR from trace
  float *yf = new float[NoBins];
  for (int i=0;i<NoBins;i++) 
    yf[i] = y[i];
  
  spec->Deconvolution(yf,SPR, NoBins, 100,20,1.5);//50,10,20.0
  
  delete [] SPR;
  delete spec;
  return yf;
};