#define ARRAY_SIZE(array) (sizeof((array))/sizeof((array[0])))

//Useful ROOT forums
//How to accurately fit a TGraph with peaks:
//https://root-forum.cern.ch/t/how-to-accurately-fit-a-tgraph-with-peaks/20063/7
//Problem with peak search tool of TSpectrum:
//https://root-forum.cern.ch/t/problem-with-peak-search-tool-of-tspectrum/13577
//TSpectrum Class Reference:
//https://root.cern.ch/doc/v608/classTSpectrum.html#a5ba181a45b117934b48c4ef5f78d0b2b
//Peaks finding:
//https://root-forum.cern.ch/t/peaks-finding/18566
//USeful for very nearby peaks:
//https://root.cern.ch/root/html534/guides/spectrum/Spectrum.html#dimensional-spectra-3

int npeaks = 10; //number of peaks to fit
double *xpeak;
double xpeakmax, a, b = 0.;
double sigmapeak[2];

//0.------Definition of the functions to be fitted to the histogram
double fpeaks(double *x, double *par) {
  double result = par[0] + par[1]*x[0];
  for (int p=0;p<npeaks;p++) {
    double norm  = par[3*p+2];
    double mean  = par[3*p+3];
    double sigma = par[3*p+4];
    result += norm*TMath::Gaus(x[0],mean,sigma);
  }
  //fpeaks = gaus_1+gaus_2+...+gaus_npeaks + linear_func.
  return result;
}

double fgauss(double *x, double *par) {
  double norm = par[0];
  double mean = par[1];
  double sigma = par[2];
  
  double y = norm*TMath::Gaus(x[0],mean,sigma);
  
  return y;
}

void findRGApeaks(){
  gStyle->SetOptStat(0);
  //gStyle->SetStatY(0.4);
  //gStyle->SetStatX(0.9);
  //gStyle->SetStatW(0.2);
  //gStyle->SetStatH(0.2);
  
  //1.--------- Fill histograms sg and bg. Subtract bg from sg.
  ifstream infile1, infile2;
  const int nmax=2560;
  float x1[nmax], x2[nmax], y1[nmax], y2[nmax];
  
  int n  = 0;
  int nn = 0;
  string line;
  TF1 *fit = new TF1();
  
  string basename = "../";
  
  string nosource = "2018.06.19_whileflushinggar_bkg_mod.txt";
  string sourcefile = "2018.06.19_whileflushinggar_mod.txt";
  
  //SIGNAL
  infile1.open(basename + sourcefile, ios::in);
  //BACKGROUND
  infile2.open(basename + nosource, ios::in);
  
  
  string outtxtfile = basename + sourcefile + "_result.txt";
  
  if (infile1.fail()){
    cout << "Sorry, couldn't open file" << endl;
    exit(1);
  }
  
  if (infile2.fail()){
    cout << "Sorry, couldn't open file" << endl;
    exit(1);
  }
  
  //Fill array, skipping empty lines
  while(getline(infile1, line)){
    infile1 >> x1[n] >> y1[n];
    n++;
  }
  
  infile1.close();
  
  while(getline(infile2, line)){
    infile2 >> x2[nn] >> y2[nn];
    nn++;
  }
  infile2.close();
  
  
  TH1F *h = new TH1F("h","sg-bg spectrum",nmax,0,80);
  TH1F *hs = new TH1F("hs","sg spectrum",nmax,0,80);
  TH1F *hb = new TH1F("hb","bg spectrum",nmax,0,80);
  
  TGraph *g1 = new TGraph(n,x1, y1);
  TGraph *g2 = new TGraph(n,x2, y2);
  
  for(int i=0; i<nmax;i++){
    if (x1[i]>27 && x1[i]<29) y1[i]*=1E-50;
    if (x2[i]>27 && x2[i]<29) y2[i]*=1E-50;
    if (x1[i]>17 && x1[i]<19) y1[i]*=1E-50;
    if (x2[i]>17 && x2[i]<19) y2[i]*=1E-50;
    //Peaks at 18 and 28 need to be reduced otherwise you get inf for that bin and hist is not drawn.
    hs->Fill(x1[i],y1[i]);
    hb->Fill(x2[i],y2[i]);
  }
  
  //cout << "size here is    " << ARRAY_SIZE(x1) << endl;
  
  h->Add(hs,hb,1,-1);

  h->Rebin(10); hs->Rebin(10); hb->Rebin(10);
  
  //2.--------- Finding peaks to fit in the spectrum
  //TSpectrum class finds candidates peaks
  TSpectrum *s = new TSpectrum(npeaks);
  int nfound = s->Search(h ,0.01,"nobackground",0.05);
  
  TPolyMarker *pm = new TPolyMarker(npeaks, s->GetPositionX(),s->GetPositionY());
  pm->SetMarkerSize(1.5);
  pm->SetMarkerColor(kRed);
  h->SetMarkerStyle(21);
  h->GetListOfFunctions()->Add(pm);
  
  printf("Found %d candidate peaks to fit\n",nfound);
  
  //Background Estimation optional
  //TF1 *fline = new TF1("fline","pol1",-2,10);// linear function (pol1) defined between [-2,10]
  //h->Fit("fline","qn");
  
  //4.-------- Setting initial parameters and fitting peaks
  if(npeaks>1) {
    double par[100]; //Array of fitting parameters
    //Loop on the found peaks for setting initial parameters. Peaks at the background level will be disregarded from the fitting.
    //fline->FixParameter(0,0);
    //fline->FixParameter(1,0);
    //par[0] = fline->GetParameter(0);
    //par[1] = fline->GetParameter(1);
    par[0] = 0;
    par[1] = 0;
    npeaks = 0;
    double *xpeaks = s->GetPositionX();//array with X-positions of the centroids found by TSpectrum
    for (int p=0;p<nfound;p++) {
      float xp = xpeaks[p];
      int bin = h->GetXaxis()->FindBin(xp);
      float yp = h->GetBinContent(bin);
      float ypb = hb->GetBinContent(bin);
      
      //Condition for valid peaks: height-sqrt(height) has to be greater than the bckgnd:
      //if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      //if (yp-TMath::Sqrt(yp) < ypb) continue;
      if (yp-TMath::Sqrt(yp) < ypb) continue;
      par[3*npeaks+2] = yp; //height
      par[3*npeaks+3] = xp; //centroid
      par[3*npeaks+4] = 0.1; //sigma -I've choosed it as same for all peaks-
      npeaks++;
    }
    printf("\n FOUND %d USEFUL PEAKS TO FIT \n",npeaks);
    cout << endl;
    //printf("Parameter No 1,2 ---> linear background Fit\n");
    printf("Parameters No 3*n+2 ---> Gauss Constants\n");
    printf("Parameters No 3*n+3 ---> Gauss Means\n");
    printf("Parameters No 3*n+4 ---> Gauss Sigmas\n");
    printf("Now fitting: Be patient\n");
    //TF1 *f = new TF1("f",fpeaks,0,10,2+3*npeaks);
    
    fit = new TF1("fit",fpeaks,0,80,3*nfound+2);
    //fit->SetNpx(10000);
    fit->SetParameters(par);
    fit->FixParameter(0,0);
    fit->FixParameter(1,0);
    
    //TVirtualFitter::Fitter(h,npeaks); //Sigma set to value above par[3*npeaks+4] = 0.1; //sigma -I've choosed it as same for all peaks-
    
    h->Fit("fit", "R"); //Sigma free to change
    
    cout << endl;
    double a = 2*TMath::Sqrt(2*TMath::Log(2));
    for(int i=0; i<npeaks; i++){
      double norm = fit->GetParameter(3*i+2);
      double normerror = fit->GetParError(3*i+2);
      double mean = fit->GetParameter(3*i+3);
      double meanerror = fit->GetParError(3*i+3);
      double sigma = fit->GetParameter(3*i+4);
      double sigmaerror = fit->GetParError(3*i+4);
      double fwhm = a*sigma;
      double fwhmerror = a*sigmaerror;
      
      cout << "Peak number " << i+1 << ": norm = " << norm << "+/-" << normerror <<
      ", mean = " << mean << "+/-" << meanerror <<
      ", sigma = " << sigma << "+/-" << sigmaerror <<
      ", FWHM = " << fwhm << "+/-" << fwhmerror << endl;
      cout << "x value as found by TSpectrum " << xpeaks[i] << endl << endl ;
    }
    
  }
  
  else {
    xpeak=s->GetPositionX();
    //cout << "xpeak : " << *(double*)xpeak << endl;
    double minx = *(double*)xpeak-0.3;
    double maxx = *(double*)xpeak+0.3;
    
    fit = new TF1("fit",fgauss,minx,maxx,3);
    fit->SetNpx(1000);
    double height = h->GetBinContent(*(double*)xpeak);
    fit->SetParameters(height, *(double*)xpeak, 0.2);
    //fit->FixParameter(1,*(double*)xpeak);
    h->Rebin(10); hs->Rebin(10); hb->Rebin(10);
    h->Fit("fit", "R");
    for (int i=0;i<3;i++) {
      printf("%g  %g\n",fit->GetParameter(i),fit->GetParError(i));
    }
  }

  //fline->SetLineColor(kGreen);
  
  h->SetLineColor(kGreen-3);
  h->SetMarkerStyle(21);
  h->SetMarkerColor(kGreen-3);
  h->SetMarkerSize(0.07);
  
  hb->SetLineColor(kBlue);
  hb->SetMarkerStyle(21);
  hb->SetMarkerColor(kBlue);
  hb->SetMarkerSize(0.07);
  
  hs->SetLineColor(kRed);
  hs->SetMarkerStyle(21);
  hs->SetMarkerColor(kRed);
  hs->SetMarkerSize(0.07);
  
  g1->SetLineColor(kRed);
  g2->SetLineColor(kBlue);
  
  TCanvas * c1 = new TCanvas("c1");
  c1->SetLogy();
  c1->cd();
  //h->Print("All");
  //hs->Draw("hist");
  h->Draw("hist");
  //hb->Draw("hist same");
  fit->Draw("same");
  h->GetYaxis()->SetTitle("Counts");
  h->GetXaxis()->SetTitle("Pulse integral [V*ns]");
  h->GetYaxis()->SetRangeUser(1E-15,1E-7);
  
  double xl1=.55, yl1=0.75, xl2=xl1+.3, yl2=yl1+.125;
  TLegend *leg = new TLegend(xl1,yl1,xl2,yl2);
  leg->AddEntry(h,"signal-bkgd","lp");
  leg->AddEntry(hs,"signal only","lp");
  leg->AddEntry(hb,"bkgd only","lp");
  leg->AddEntry(fit,"sum of gauss fit","lp");
  leg->Draw("same");
  
  /*
   TCanvas * c2 = new TCanvas("c2");
  c2->SetLogy();
  c2->cd();
  g1->Draw();
  g2->Draw("same");
  */
  
}