#include <math.h>

TH1F *autocorr(const TH1F *hist)
{
   Int_t nbins = hist->GetNbinsX();
   std::cout << "number of bins is " << nbins << std::endl;
   Double_t *seq = new Double_t[nbins];
   for (int bin=0; bin < nbins; ++bin) {
      seq[bin] = hist->GetBinContent(bin);
   }
      
   TH1F *hcor = new TH1F("acor","autocorrelation plot",1000,0,1000);

   for (int n=0; n < 1000; ++n) {
      Double_t sum=0;
      for (int bin=0; bin+n < nbins; ++bin) {
         sum += seq[bin]*seq[bin+n];
      }
      hcor->Fill(n+0.5,sum);
   }

   delete [] seq;
   return hcor;
}

TH1F *foldit(const TH1F *hist)
{
   Double_t period = 35.78;
   Int_t nbins = hist->GetNbinsX();
   std::cout << "number of input bins is " << nbins << std::endl;
   Double_t xmin = hist->GetXaxis()->GetXmin();
   Double_t xmax = hist->GetXaxis()->GetXmax();

   Int_t nbout = period*10+0.99;
   std::cout << "number of output bins is " << nbout << std::endl;
   Double_t xrange = (xmax-xmin)*nbout/nbins;
   TH1F *hout = new TH1F("acor","folded t-tRF difference (ns)",
                         nbout,0,xrange);
   for (int bin=0; bin < nbins; ++bin) {
      Double_t x = ((bin % nbout) + 0.5)*xrange/nbout;
      hout->Fill(x,hist->GetBinContent(bin));
   }
   return hout;
}

TF1 *fitf(Double_t period=2.05,
          Double_t height=50,
          Double_t sigma=0.35,
          Double_t origin=-0.8)
{
   TF1 *f = new TF1("fitf", fitfun, 0, 25, 4);
   f->SetParameters(period,height,sigma,origin);
   f->SetParNames("period","height","sigma","origin");
   return f;
}

double sqr(double x) { return x*x; }

Double_t fitfun(Double_t *x, Double_t *params)
{
   Double_t period = params[0];
   Double_t height = params[1];
   Double_t sigma = params[2];
   Double_t origin = params[3];

   Double_t sum = 0;
   for (int peak=-5; peak < 25; ++peak) {
     sum += height*exp(-0.5*sqr((x[0]-period*peak-origin)/sigma));
   }
   return sum;
}