#include "TSpectrum"
#include "TSpectrumFit"
#include "TH1.h"
#include "TF1.h"
#include "TVirtualFitter.h"

#include <stdio.h>
#include <iostream>
#include <fstream>

Int_t npeaks = 30;
Float_t factor=1.e-6*1.e-26/1e4;//factor of 10^-6 due to micro Jy, factor of 1e-26 to conversion to J s-1 cm-2 Hz-1// in order to have cm^-2 its 10^4 cm^2
Float_t c=2.998e8; //speed of light in vaccuum in m/s reftractive index of air would be 1.000277
Float_t planck=6.62606957e-34; //planck constant to calculate enrgy of a photon with given frequency
Float_t FoV=12600.;//MAGIC FoV of 3.5 deg expressed in arsec
Float_t area=236e4;//Area of one MAGIC telescope in cm^2
//generate n peaks at random
   Double_t par[3000];
   par[0] = 0.8;
   par[1] = -0.6/1000;
Double_t fpeaks(Double_t *x, Double_t *par) {
   Double_t result = par[0] + par[1]*x[0];
   for (Int_t p=0;p<npeaks;p++) {
      Double_t norm  = par[3*p+2];
      Double_t mean  = par[3*p+3];
      Double_t sigma = par[3*p+4];
      result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
}

void FitPeaks(TString file,Int_t np=40){
  gStyle->SetOptStat(0);
  Double_t x,y;
 TNtupleD* data = new TNtupleD("data", "data", "x:y");
  data->ReadFile(file);
  // if (file==NULL){
  // printf("File %s not found\n",data);
  // break;
  // }
TCanvas *ntup=new TCanvas("ntuple", "ntuple");
  cout<<"Read "<<file<<endl;
  data->SetBranchAddress("x",&x);
  data->SetBranchAddress("y",&y);
   data->Draw("y:x","","l");
   TF1 *fbg=new TF1("fbg","1.42766e-2*x",400,1000);
   
  Int_t  n = data->GetEntries();
  Float_t nu;
  Double_t e;
  cout<<"Spectrum contains "<<n<<" points"<<endl;
  data->GetEntry(0);
  Float_t xmin=x;
  data->GetEntry(n-1);
  Float_t xmax=x;
  
  TGraph *g=(TGraph*)ntup->GetListOfPrimitives()->FindObject("Graph");
  
  

 TH2F *temp=(TH2F*)ntup->GetListOfPrimitives()->FindObject("htemp");
 TH1F *h=new TH1F();
  TH1F *h1= new TH1F("h1","h1",n,htemp->GetXaxis()->GetXmin(),htemp->GetXaxis()->GetXmax());
  TH1F *h2= new TH1F("h2","h2",n,htemp->GetXaxis()->GetXmin(),htemp->GetXaxis()->GetXmax());


 
  for(Int_t i=0;i<n;i++){
    data->GetEvent(i);
    nu=c/(x*1.e-9);
    e=nu*planck;
    h->Fill(x,y);	
    h1->Fill(x,y/e*factor);
    h2->Fill(x,y/e*factor*area*FoV**2);
  }
  TCanvas *c4=new TCanvas();
 h->Draw("same l");
break;
   // c1->Update();
  TCanvas *c=new TCanvas("c","c",1000,1000);
  c->Divide(1,4);
  c->cd(1);
  h->GetYaxis()->SetTitle("Intensity [#muJy arcsec^{-2}] ");
  h->GetXaxis()->SetTitle("#lambda [nm]");
  h->Draw();
     c->cd(2);
    h1->GetYaxis()->SetTitle("Intensity [photons s^{-1} cm^{-2} arcsec^{-2}]");
    h1->GetXaxis()->SetTitle("#lambda [nm]");
    h1->Draw();

    
   
  
     c->cd(3);
    h2->GetYaxis()->SetTitle("Intensity [photons s^{-1}]");
    h2->GetXaxis()->SetTitle("#lambda [nm]");
      h2->Draw();
      c->cd(4);
  npeaks = TMath::Abs(np);
  TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
  f->SetNpx(1000);
   TH1F *h3 = (TH1F*)h->Clone("h3");
   //Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h,1,"",0.000001);
   printf("Found %d candidate peaks to fit\n",nfound);
   //Estimate background using TSpectrum::Background
   TH1 *hb = s->Background(h,30,"same");
   if (hb) c->Update();
   if (np <0) return;
 //estimate linear background using a fitting method
 
   TF1 *fline = new TF1("fline","pol1",0,1000);
   h3->Fit("fline","qn");
   //Loop on all found peaks. Eliminate peaks at the background level
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);
   npeaks = 0;
   Float_t *xpeaks = s->GetPositionX();
   for (Int_t p=0;p<nfound;p++) {
      Float_t xp = xpeaks[p];
      Int_t bin = h3->GetXaxis()->FindBin(xp);
      Float_t yp = h3->GetBinContent(bin);
      if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      par[3*npeaks+2] = yp;
      par[3*npeaks+3] = xp;
      par[3*npeaks+4] = 3;
      npeaks++;
   }
   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Now fitting: Be patient\n");
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
   //we may have more than the default 25 parameters
   TVirtualFitter::Fitter(h3,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   h3->Fit("fit");  
   c->cd();
   //c->SaveAs("NSB_Fit.root");
   //c->SaveAs("NSB_Fit.C");
 
}